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1.  INTRODUCTION 


The  DIRTRAN-I  Code  Is  a  computer-implemented  model  for  predicting 
the  optical  effects  of  an  explosion-produced  dust  cloud  as  it  disperses  in 
the  lower  atmosphere.  This  model  is  based  on  first  principles  of  fluid 
dynamics,  atmospherics,  and  optics.  The  model  has  been  validated  using 
cloud  dimension  and  line-of-sight  optical  transmission  data  from  the 
DIRT  I  and  Graf  II-Winter  Army  dust  obscuration  field  trials.  The 
DIRTRAN-I  Code  exploits  information  available  about  crater  sizes  produced 
by  explosions  in  conjunction  with  distinct  models  for  coupling  of  energy 
to  the  ground  for  artillery  projectiles  versus  bare  charges.  The  model 
recognizes  that  dust  ejecta  are  partitioned  into  a  buoyantly  rising  fire¬ 
ball  and  a  non-buoyant  "dust  skirt"  which  is  subject  to  diffusion  in  the 
vertically  sheared  wind  field.  The  DIRTRAN-I  Code  solves  separately  for 
these  two  clouds.  The  solutions  are  based  on  atmospheric  diffusion  theory 
and  take  into  account  the  effects  of  wind  and  temperature  profiles  in  the 
constant  shear  stress  layer  of  the  lower  atmosphere  for  different  atmo¬ 
spheric  stability  categories.  Separate  treatment  is  given  to  particles  of 
different  sizes,  the  larger  ones  being  allowed  to  settle  out.  Outputs  of 
the  code  include  dust  cloud  displacement  and  dimensions  for  both  the  non- 
buoyant  wind-dominated  skirt  and  the  initial  buoyant  fireball  as  it  is  wind 
blown  and  eventually  also  becomes  subject  to  wind  diffusion.  Line-of-sight 
transmittances  at  several  wavelength  bands  (visible:  0.4  -  0.7  pm; 
infrared:  0.8  -  1.1,  3.5  -  4.0,  and  8.5  -  12  pm;  mm  wave:  94  -  140  GHz) 
are  also  output  options. 

Because  of  the  analytic  solutions  derived  from  first-principle  physics, 
the  DIRTRAN-I  Code  has  been  designed  to  keep  both  storage  and  computation 
time  to  a  minimum.  In  order  to  be  machine  transportable,  the  code  has 
been  written  in  ANSI  FORTRAN  IV  with  no  system  specific  enhancements. 
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2.  FUNCTIONS  AMD  SUBROUTINES 


There 
Section  5. 
ture  of  the 


is  a  great  deal  of  internal  documentation 
In  Figs.  2.1  through  2.5  are  box  diagrams 
code . 


in  the  code  listed  in 
illustrating  the  struc- 


Figure  2.1  Controlling  Routines. 
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Figure  2.4  Routines  for  Wind  Dispersion  Model. 


COMPUTES  VISIBLE  CLOUD  DIMENSIONS 
FROM  BUOYANT  AND  WIND  DISPERSED 


Figure  2.5  Routines  for  Determining  Cloud  Dimensions. 


Miscellaneous 

FUNCT 

GRAD2 

GFUN 

PERP 

UNIT 

VS  UN 


Routines 

Two-dimensional  function  which  provides  optically  weighted 
column  densities  for  visible  contour  tracing. 

Computes  the  two-dimensional  gradient  vector  of  a  function. 

The  restriction  of  a  function  of  two  variables  to  a  line 
specified  by  a  base  point  and  direction  vector. 

Computes  a  vector  rotated  90°  counterclockwise  from  a 
given  vector. 

Computes  the  length  of  a  vector  and  a  unit  vector  in  the 
same  direction  as  a  given  vector. 

Adds  a  given  scaler  multiple  of  one  given  vector  to  another. 
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3.  OSAGE 


The  user  accesses  the  DIRTRAN-I  Code  via  Subroutine  DRTRAN. 

CALL  DRTRAN  (WAVE1,  ICLMAT,  TRNLOS ,  IERR) 

where : 

WAVE1  is  the  wavelength  in  micrometers.  Valid  ranges 

are:  0.4  -  0.7  micrometers 

0.8  -  1.1  micrometers 

3.5  -  4.0  micrometers 

8.5  -  12.0  micrometers 

2100.0  -  3200.0  micrometers 

ICLMAT  is  an  integer  index  which  is  0  if  meteorological 
data  is  read  from  an  external  unit  with  other 
inputs  and  is  1  if  passed  in  COMMON  /CLYMAT/. 

TRNLOS  is  the  returned  value  of  the  transmittance  along 
the  line-of-sight  between  the  transmitter  and  the 
rec iever . 

IERR  is  an’  integer  error  code  which  is  returned  one  if 
a  fatal  error  occurs  and  zero  otherwise. 

3. 1  Common  Blocks  Used  to  Pass  Information  to  DIRTRAN-I 

There  are  two  common  blocks  which  communicate  information  to  DIRTRAN-I: 

/I0UNIT/  and  /CLYMAT/. 

3.1.1  COMMON  /IOUNIT/  IOIN.  I00UT.  ISPTPF.  LOUNIT.  NDIRTU.  NBSCAT 

Of  the  variables  in  /IOUNIT/  only  IOIN,  IOOUT,  and  NDIRTU  are  used  by  DIRTRAN-1. 

IOIN  The  Fortran  logical  unit  from  which  DIRTRAN-I 

reads  input  data  to  be  described  in  Subsection  3.2. 
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IOOUT  The  Fortran  logical  unit  onto  which  DIRTRAN-I 
writes  output  and  error  messages. 

NDIRTU  The  Fortran  logical  unit  from  which  DIRTRAN-I 

reads  the  data  file  containing  tabulated  values 
for  the  moments  of  the  atmospheric  diffusion 
equation . 


3.1.2  COMMON  /CLYMAT/  TEMP.  PRESS.  RH.  AH.  DP.  VIS.  CLDAMT.  CLDHYT.  FOGPRB . 


IPASCT .  WNDVEL .  WNDDIR 


DIRTRAN-I  uses  only  TEMP,  IPASCT,  WNDVEL,  and  WNDDIR. 


TEMP  Temperature  in  degrees  C  taken  at  approximately 
two  meters  above  ground. 


IPASCT  Integer  1-6  corresponding  to  Pasquill 
Categories  A  -  F. 


WNDVEL  The  wind  speed  in  meters  per  second  measured 
at  approximately  two  meters  above  ground. 

WNDDIR  The  wind  direction  (  in  degrees  clockwise  from 
true  north)  from  which  the  wind  is  blowing. 

With  this  option  the  user's  coordinate  system 
must  have  the  positive  x-axis  pointing  east  and 
positive  y-axis  pointing  north.  (Thus  0°  corresponds 
to  a  wind  blowing  from  the  north  to  the  south;  90°  is 
a  wind  blowing  from  the  east  to  the  west;  etc.) 


Input  from  FORTRAN  Logical  Unit  JOIN 


The  input  read  from  logical  unit,  IOIN,  provides  the  user  x^ith  a  variety 
of  ways  to  use  DIRTRAN-1.  There  are  seven  types  of  records  to  be  distinguished 
on  the  input  file.  They  are  listed  here  with  the  names  of  variables  contained 
on  each  record  and  the  format  type. 
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RECORD 


NAME 

NEWATM 

NEWSRC 

LOSTRN 

EDGE 

NEWTIM 

FORMAT 

LI 

LI 

LI 

LI 

LI 

NAME 

NATMOS 

ZTMPj 

TMPMES 

ZWND 

WNDMES 

ZTMP 

TEMPMESl 

ZWND 

WNDMES 

ZINV 

THWND 

FORMAT 

11, IX 

CM 

r>* 

Pn 

F7.2 

F7.2 

F7.2 

F7.2 

F7.2  1 

F7 . 2 

F7.2 

F7.2 

F7.2 

NAME  NSOIL  NCHRG  CHWT 

IDETDEP  DSOD 

ISRCCOR  SRCCOR 

FORMAT  11, IX  11, IX  F7.2 

F7.2  F7.2 

f  F7.2  F7.2 

NAME 

TRNCOR 

TRNCOR 

TRNCOR 

RECCOR 

RECCOR 

RECCOR 

FORMAT 

F7.2 

F7.2 

F7.2 

F7.2 

F7.2 

F7.2 

NAME 

lOBSCOR 

OBSCOR 

SPCHT 

FORMAT 

I  F7.2 

F7.2 

F7.2 

NAME 

FORMAT 

Itime 

IF7.2 

NAME 

FORMAT 

ICNTNU 

1  LI 

The  description  of  the  variables  is  in  the  comments  in  subroutine  DRTRAN  and 
appears  on  Pages  5-1  -  5-6  in  this  manual.  Additional  descriptions  for  some 
of  the  variables  appears  in  Section  4. 

The  input  file  contains  one  or  several  sequences  of  these  records  each 
of  which  must  begin  with  record  1,  end  with  record  7,  and  contain  a  subset  of 
records  2-6  corresponding  to  the  entries  in  record  1.  In  each  sequence, 
each  of  records  2  through  6  must  appear  if  and  only  if  the  corresponding 
control  variable  entered  in  record  1  is  .TRUE.. 

If  NEWATM  is  .TRUE,  then  record  2  must  appear. 

If  NEWSRC  is  .TRUE,  then  record  3  must  appear. 

If  LOSTRN  is  .TRUE,  then  record  4  must  appear. 

If  EDGE  is  .TRUE,  then  record  5  must  appear. 

If  NEWTIM  is  .TRUE,  then  record  6  must  appear. 

On  the  first  record  of  the  input  file  NEWATM,  NEWSRC,  and  NEWTIM  must 
all  be  .TRUE..  This  initializes  the  meteorological  conditions,  the  charge  3nd 
soil  characteristics,  and  time  after  blast  for  the  first  observation.  After 
that,  DIRTRAN-I  assumes  that: 


i 


i 


1 

I 
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•  Meteorological  data  are  unchanged  until  NEWATM  is  .TRUE, 
again. 

•  Charge  and  Soil  Characteristics  are  unchanged  until 
NEWSRC  is  .TRUE,  again. 

•  Time  after  blast  is  unchanged  until  NEWTIM  is  .TRUE,  again. 

The  first  time  that  NEWATM  is  .TRUE.,  if  ICLMAT  ■  1  ( see  arguments  of 
DRTRAN  in  Section  2),  then  record  2  should  not  appear  following  record  1. 
This  is  the  only  time  that  a  control  variable  in  record  1  may  be  .TRUE, 
without  the  corresponding  record  2-6  following. 

Record  7  contains  the  control  variable,  CNTNU,  which  is  .TRUE,  if 
another  sequence  of  records  1  -  7  is  to  be  read  on  this  call  of  DRTRAN  and 
.FALSE,  is  DRTRAN  is  to  return  control  to  the  calling  routine.  An  example 
follows . 

3 . 3  Example 


T 

FFTTT 
175.00- 
175. 0- 
2. 

T 

FFTTT 
1 75.  DO- 
175.  0- 
5. 

T 

FFTTT 
175. 00- 
175.0- 
10. 

T 

FFTTT 
1 75. 00- 
175.0- 
20. 

T 

FFTTT 
175.00- 
175.  0- 
30. 

T 

FFTTT 

175.00- 

175.0- 

40. 

T 

FFTTT 

175.00- 

175.0- 

60. 

F 


1148.9 

1148.0 


1148.9 

1148.0 


1148.9 

1148.0 


1148.9 

1148.0 


1148.9 

1148.0 


1148.9 

1148.0 


1148.9 

1148.0 


1.8  -173.0  1135.0 

1.8 


1.8  -173.0  1135.0 

1.8 


1.8  -173.0  1135.0 

1.8 


1.8  -173.0  1135.0 

1.8 


1.8  -173.0  1135.0 
1*8 


1.8  -173.0  1135.0 

1.8 


1.8  -173.0  1135.0 
1.  8 


1.8 


1.8 


1.8 


1.8 


1.8 


1.8 


1.8 


The  results  from  the  above  input  file  follow 


DIRTRAN-1  DUST  CLOUD  INFRARED  TRANSMISSION  CALCULATION 


ALL  UNITS  ARE  MRS  UNLESS  OTHERWISE  SPECIFIED. 


PASQU ILL  CATEGORY 

HT  1.00  TEMP 
WIND  DIRECTION 

SOIL  INDEX  1 

CHARGE  INDEX  1 
WEI  CUT  OF  CHARGE 
DETONATION  DEPTH 
DEPTH  OF  SOD 
SOURCE  COORDINATES 


3 

2  71.30  HT 
90.  00 


15.00 
0.  00 
0.  1  5 
0. 


1.0  0  WIND  2.3  0 


0.  00 


TIME  AFTER  BLAST 


1.  OC 


WAVELENGTH 

TRANSMITTER  COORDINATES 
RECEIVER  COORDINATES 
TRANSMITTANCE  ALONG  THE  LINE 


0. 55  MICROMETERS 
175.00  -1148.00 

-1  73.  00  1  1  35.  00 

-OF-SIGHT  0.149E-02 


1.  80 

1.  PO 


OBSERVER  COORDINATES  175.00 

THE  HEIGHT  OF  THE  CLOUD  IS 
THE  CENTROID  COORDINATES  ARE 
THE  WIDTH  AT  THE  CENTROID  IS 

THE  WIDTH  AT  1.80  METERS  IS 

5  CONTOUR  POINTS  HAVE  BEEN  DETERMINED 
-6.654  1.800 

-3.1G4  4.794 

0.410  6. 388 

4.  003  4.  794 

7.400  1.800 


I  148. 00 
8.39  METERS 
0.41  4.79 

7.  1  9  METEPS 
14.06  METERS 


TIME  AFTER  BLAST 


2.  00 


WAVELENGTH 

TRANSMITTER  COORDINATES 
RECEIVER  COORDINATES 
TRANSMITTANCE  ALONG  THE 


0.55  MICROMETERS 
175.00  -1148.90 

-1  73.  00  1  1  35.  00 

LINE-OF-S ICHT  0.296E-02 


1.80 
1. 80 


OBSERVER  COORDINATES  175.00  -1148.00 

THE  HEIGHT  OF  THE  CLOUD  IS  12.83  METERS 

THE  CENTROID  COORDINATES  ARE  0.83  e.34 
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THE  WIDTH  AT  THE  CENTROID  IS  8.97  METERS 

THE  WIDTH  AT  1.80  METERS  IS  14.69  METERS 

5  CONTOUR  POINTS  HAVE  BEEN  DETERMINED 
-6.590  1.800 

-3.650  8.345 

0.832  12.827 

5.  315  8.  345 

8.  098  1.800 


TIME  AFTER  BLAST 
WAVELENGTH 

TRANSMITTER  COORDINATES 
RECEIVER  COORDINATES 
TRANSMITTANCE  ALONG  THE  LINE 

OBSERVER  COORDINATES 
THE  HEIC1IT  OF  THE  CLOUD  IS 
THE  CENTROID  COORDINATES  ARE 
THE  WIDTH.  AT  THE  CENTROID  IS 
THE  WIDTH  AT  1.80  METERS 

5  CONTOUR  POINTS  HAVE  BEEN 
-6.396  1.800 

-4.409  16.845 

2.  206  23.  462 

£.825  16.845 

10.166  1.800 


5.00 

0.55  M  ICROMETERS 
1  75.  00  -1  148.  90  1.80 

-1  73.  00  1  1  35.00  1.80 

OF-SICHT  0.156E-01 

175.00  -1148.00 

23.  46  METERS 
2.21  16.85 

1  3.  2  3  METERS 
IS  16.56  METERS 

DETERMINED 


TIME  AFTER  BLAST  10.00 

WAVELENGTH  0. 5 5  M I CROM ETER S 

TRANSMITTER  COORDINATES  175.00  -1148.90  1.80 

RECEIVER  COORDINATES  -173.00  1135.00  1.80 

TRANSMITTANCE  ALONG  THE  L IN E-0 F-S ICHT  0.I07E  00 


OBSERVER  COORDINATES  175.00 

THE  HEIGHT  OF  THE  CLOUD  IS 
THE  CENTROID  COORDINATES  ARE 
THE  WIDTH  AT  THE  CENTROID  IS 
THE  WIDTH  AT  1.80  METERS  IS 

5  CONTOUR  POINTS  HAVE  BEEN  DETERMINED 
-5.  762  1.800 

-4.453  26.943 

4.  709  36.  1  05 

13.871  26.943 

13.613  1.800 


-1  148.00 
36.  1  0  METERS 
4.71  26.94 

1 6. 32  ME7EPS 
1  9.  37  METERS 
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TIME  AFTER  BLAST 


20.  00 


WAVELENGTH 

TRANSMITTER  COORDINATES 
RECEIVER  COORDINATES 
TRANSMITTANCE  ALONG  THE 


0. 55  MICROMETERS 
175.00  -1148.90 

-173.00  1135.00 

LINE-OF-S ICHT  0. 508E  00 


1.80 

1.80 


OBSERVER  COORDINATES 
THE  HEIGHT  OF  THE  CLOUD  IS 
THE  CENTROID  COORDINATES  ARE 
THE  WIDTH  AT  THE  CENTROID  IS 
THE  WIDTH  AT  1.80  METERS 

5  CONTOUR  POINTS  HAVE  BEEN 


175. 00 


IS 

DETERMINED 


- 

4. 

1 

79 

1. 

800 

-1 

1. 

0 

30 

40. 

057 

1 

0. 

0 

32 

59. 

232 

3 

0. 

9 

83 

40. 

057 

1 

9. 

8 

S3 

1. 

800 

-1148.00 
59.  2  3  METERS 
9.  98  40. 

42.  01  KETEPS 
24.  06  METERS 


06 


TIME  AFTER  BLAST 


30.  00 


WAVELENGTH 

TRANSMITTER  COORDINATES 
RECEIVER  COORDINATES 
TRANSMITTANCE  ALONC  THE 


0.55  MICROMETERS 
175.00  -1148.90 

-1  73.  00  1  1  35.  00 

LINE-OF-S  ICHT  0.  789E  00 


1.  FO 
1.  FO 


OBSERVER  COORDINATES 
THE  HEIGHT  OF  THE  CLOUD  IS 
THE  CENTROID  COORDINATES  ARE 
THE  WIDTH  AT  THE  CENTROID  IS 
THE  WIDTH  AT  1.80  METERS 

5  CONTOUR  POINTS  HAVE  BEEN 
-2.285  1.E0C 

-16.638  40.178 

1  7.  039  67  .  7  79 

50.  4  1  3  40.  1  78 

25.528  1.800 


1  75.  00 


IS 

DETERMINED 


■1  146.00 
6  7.  7  8  METERS 

1  5.  89  40.  1  8 

69.  05  METERS 

2  7.81  METERS 


TIME  AFTER  BLAST 


40.  00 


WAVELENGTH 

TRANSMITTER  COORDINATES 
RECEIVER  COORDINATES 
TRANSMITTANCE  ALONG  THE 


0.  55  MICROMETERS 
1  75.  00  -1  148.  90 

-173.00  1135.00 

LINE-OF-S ICHT  0. 909E  00 


i.  eo 

1.  FO 


1  t 


OBSERVER  COORDINATES  175.00 

THE  HEIGHT  OF  THE  CLOUD  IS 
THE  CENTROID  COORDINATES  ARE 
THE  WIDTH  AT  THE  CENTROID  IS 
THE  WIDTH  AT  1.80  METERS  IS 

5  CONTOUR  POINTS  HAVE  BEEN  DETERMINED 
0.  235  1.800 

-19.906  40.703 

20.352  72.355 

63.238  40.703 

30.860  1.800 


-1 148.00 
72. 35  METERS 
2  1.  6  7  40.  70 

83.  14  METERS 
30.  63  METERS 


TIME  AFTER  BLAST  60.00 

WAVELENGTH  0.55  MICROMETERS 

TRANSMITTER  COORDINATES  175.00  -1148.90  1.80 

RECEIVER  COORDINATES  -173.00  1135.00  1.80 

TRANSMITTANCE  ALONG  THE  L I N E -0 F-S IGHT  0.977E  00 

OBSERVER  COORDINATES  175.00  -1148. CO 

THE  HEIGHT  OF  THE  CLOUD  IS  77.53  METERS 

THE  CENTROID  COORDINATES  ARE  33.99  46.28 

THE  WIDTH  AT  THE  CENTROID  IS  99.13  METERS 

THE  WIDTH  AT  1.80  METERS  IS  36.56  METERS 

5  CONTOUR  POINTS  HAVE  BEEN  DETERMINED 
4.962  1.800 

-1  5.  5  7  1  48.  2  84 

34.  569  7  7.  5  35 

83.556  48.284 

41.524  1 . 800 
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H.  DIRTRAN-I  GEOMETRY 


4.1  User  Coordinates 

A  top  view  of  the  (x,y)  user  coordinate  plane  for  the  DIRTRAN-I  geometry 
is  shown  in  Fig.  4.1.  All  inputs  are  in  this  user  coordinate  plane.  Negative 
coordinates  are  allowed. 

4.2  Wind  Geometry 

THWND  is  the  angle  in  degrees  that  the  wind  velocity  vector  makes  with 
the  user's  positive  x-axis  and  is  measured  counter-clockwise.  If  the  x-axis 
points  east,  then  the  following  relationship  holds  between  THWND  and  WINDIR 
(which  is  the  standard  meteorological  wind  direction  as  discussed  in  Sub- 
Section  3.1.2)  : 

THWND  -  270  -  WINDIR 


4.3  Transmission  Geometry 

The  height  above  ground  of  the  transmitter  and  receiver,  TRNCOR(3)  and 
RECC0R(3),  must  be  equal  and  must  be  between  0  and  5  meters.  The  transmitter 
height  will  be  used  as  a  default  for  the  receiver  height. 

4 . 4  Observer  Geometry 

A  local  coordinate  system,  x',  y',  z'  is  used  for  viewing  geometry. 

The  x'-axis  is  aligned  from  the  detonation  point  with  coordinates,  SRCCOR(I), 
to  the  observer  position  with  coordinates,  OBSCOR(I).  The  y'-axis  is  rotated 
90°  counterclockwise  from  the  x'-axis, and  the  z'-axis  is  vertical  so  as  to  make 
a  right-handed  coordinate  system.  The  visible  cloud  is  projected  along  the 
observer’s  line-of-sight,  the  x’-axis.  The  visible  cloud  contour  thus 
appears  on  the  y'-z'  coordinate  plane,  shown  in  Fig.  4.2.  Five  contour  points 
are  reported  by  DIRTRAN-I: 


t  i 
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Figure  4.2  Observer’s  View. 


•  The  left  and  right  edges  of  the  cloud  at  a  specified  height,  SPCHT, 

•  The  left  and  right  edges  of  the  cloud  at  the  centroid  height,  and 

•  The  top  of  the  cloud. 


The  centroid  is  determined  by  first  locating  the  leading  edge  of  the 
cloud  which  is  (CPTS(1,4),  CPTS(2,4))  in  Fig.  4.2.  A  horizontal  line  is 
drawn  across  the  cloud  to  determine  (CPTS(1,2),  CPTS(2,2))  in  Fig.  4.2.  The 
centroid  is  defined  to  be  the  midpoint  of  these  two. 
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5.  LISTING  OF  DIRTRAN-I  CODE 


C  DIRTRAN-I  CODE 

C 

C  TEST  CALLING  ROUTINE 

C 

COMMON  /IOUNIT/  1 OIN , I OOUT , I SPTPF , LOUNIT , NDIRTU , NB SC AT 
DATA  I01N, IOOUT.NDIRTU/5, 1, 7/ 

WAVE1-.  5  5 
ICLMAT-0 

CALL  DRTRAN (WAVE  1, ICLMAT, TRNLOS , I  ERR ) 

STOP 

END 

C  UTILITY  ROUTINE  FOR  INTERFACING  DIRTRAN-I  WITH  EC-SAEL 

C 

C  *******************  SUBFILE  1  ***************************** 

C 

SUBROUTINE  DRTRAN (WAVE  1 , ICLMAT, TRNLOS, I  ERR) 

IMPLICIT  INTEGER**  (I-N) 

LOGICAL  NEWATM,NEWSRC,LOSTRN, EDGE , NEWT IM , CNTNU, CLMRED 
DIMENSION  ZTMP(2),TMPMES(2) , ZWND ( 2 ) , WNDM ES ( 2 ) , SRCCOR ( 2 ) ,TRNCOR(3) 
1  , RECCOR (3 ) ,CPTS(2, 200),CNTRD(2),OBSCOR(2) 

COMMON  /IOUNIT/  I  0 1 N , I OOUT , I SPTPF .LOUNIT , NDIRTU , NBSC AT 
COMMON  /CLYMAT/  TEM P, PRESS , RH, AH, D P, V I S , CLDAMT , CLDH YT , 

1  FOGPRB.IPASCT.WNDVEL.WINDIR 

C  ******************************************************************** 

C 

C  DAVID  DVORE 

C  AERODYNE  RESEARCH,  INC. 

C  (6  17  )2  75-9400  X  1  27 

C 

C  PURPOSE 

C 

C  DRTRAN  CALCULATES  DUST  CLOUD  DIMENSIONS  AND  TRANSM I TTANC  ES 

C  THROUGH  DUST  CLOUDS  FOR  GIVEN  METEOROLOGICAL  DATA,  SOIL  TYPE, 

C  BOMB  CHARACTERISTICS,  AND  WAVELENCTH. 

C 

C 

C  INPUT 

NEW.NTM  A  LOGICAL  VARIABLE  WHICH  IS  .TRUE.  IF  NEW  ATMOSPHERIC 

CONDITIONS  ARE  TO  BE  SPECIFIED  AND  .FALSE.  IF  PREVIOUS 
VALUES  ARE  TO  BE  ASSUMED.  IF  .TRUE.  THEN  NATMOS,  ZTMP, 
TMPMES,  ZWND,  WNDMES,  HTINV,  AND  THWND  MUST  BE 
SPECIFIED 
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NATMOS 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ZTMP 

TMPMES 

ZWND 

WNDMES 

HTIN  V 

THWND 

NEW SRC 

CHWT 

NCHRG 


INTEGER  WITH  VALUES  1  TO  6  CORRESPONDING  TO  PASQUILL 
CATEGORIES  A  TO  F.  IF  WIND  AND  TEMPERATURE 
MEASUREMENTS  ARE  AVAILABLE  AT  TWO  HEIGHTS  THEN  NATMOS 
SHOULD  BE  SPECIFIED  0  . 

A  SINGLY  DIMENSIONED  ARRAY  CONTAINING  TWO  HEIGHTS, 
MEASURED  IN  METERS,  AT  WHICH  TEMPERATURE  MEASUREMENTS 
ARE  AVAILABLE.  THE  HEIGHTS  SHOULD  BE  IN  INCREASING 
ORDER.  IF  ONLY  ONE  IS  AVAILABLE,  THE  PASQUILL  CATEGORY 
MUST  BE  SPECIFIED  IN  NATMOS.  VALID  RANGE:  0.5-100.  M. 

A  SINGLY  DIMENSIONED  ARRAY  CONTAINING  ONE  OR  TWO 
TEMPERATURE  MEASUREMENTS  IN  DEGREES  KELVIN  TAKEN  AT 
HEIGHTS  ZTMP.  VALID  RANGE:  2  73.-31  5.  K. 

A  SINGLY  DIMENSIONED  ARRAY  CONTAINING  TWO  HEIGHTS, 
MEASURED  IN  METERS,  AT  WHICH  WIND  MEASUREMENTS  ARE 
AVAILABLE.  THE  HEIGHTS  SHOULD  BE  IN  INCREASING  ORDER. 
IF  ONLY  ONE  IS  AVAILABLE,  THE  PASQUILL  CATEGORY  MUST 
BE  SPECIFIED  IN  NATMOS.  VALID  RANCE:  0. 5-100. M. 

A  SINGLY  DIMENSIONED  ARRAY  CONTAINING  WIND  SPEEDS 
MEASURED  IN  METERS/SEC,  AT  HEIGHTS  ZWND. 

VALID  RANGE:  0.1  -  20.  M/S. 

THE  INVERSION  HEIGHT  MEASURED  IN  METERS.  DEFAULTS 
TO  200  METERS  IF  A  NUMBER  LESS  THAN  10.0  IS  GIVEN. 
VALID  RANGE:  10.0  -  200.0  M. 

THE  ANCLE  THAT  THE  WIND  VELOCITY  VECTOR  MAKES 
WITH  THE  POSITIVE  X  AXIS  MEASURED  IN  DECREES 
COUNTERCLOCKWISE. 

VALID  RANGE:  -360.  0  -  360.  0  DEGREES. 

A  LOGICAL  VARIBLE  WHICH  IS  .TRUE.  IF  A  NEW  SOURCE 
IS  TO  BE  SPECIFIED  AND  .FALSE.  IF  FURTHER  RESULTS 
ARE  DESIRED  FROM  THE  PRESENT  SOURCE.  IF  .TRUE. 

THEN  CHWT,  NCHRG,  DETDEP,  NSOIL,  DSOD,  NWL ,  AND  SRCCOR 
MUST  BE  SPECIFIED 

THE  WEIGHT  OF  THE  CHARGE  IN  LBS-TNT. 

VALID  RANGE:  1.0  -  150.0  LBS-TNT. 


CHARGE  TYPE  INDEX  WITH  FOLLOWING  VALUES 

1  SURFACE  -  LIVE  FIRE  OR  30  DECREE  TILTED 
STATIC,  TIP  ON  GROUND 

2  BARE  CHARGE  ON  SURFACE 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


3  30  DEGREE  TILTED  TIP  AT  0.3  METER  DEPTH 

4  30  DEGREE  TILTED  TIP  AT  0.6  METER  DEPTH 

5  HORIZONTAL  PROJECTILE  ON  SURFACE 

DEFAULT  VALUE  IS  1  IF  NCHRG  IS  NOT  BETWEEN  1  AND  5. 


DETDEP  THE  DEPTH  OF  DETONATION  IN  METERS. 

VALID  RANGE:  0.0  -  2.0  M. 

NSOIL  INTEGER  INDEX  OF  SOIL  TYPE.  NSOIL  IS 

1  FOR  EUROPEAN  SOIL 

2  FOR  DESERT  SOIL 

DSOD  DEPTH  OF  SOD  IN  METERS 

VALID  RANGE:  0.0  -  1.0  M. 

WAVE  1  WAVELENGTH  IN  MICROMETERS.  USED  TO  DETERMINE  NWL . 

VALID  RANGES: 


INTEGER 

INDEX 

FOR 

WAVELENGTH 

1 

FOR 

0.  4 

-0.7  M  ICROMETER 

2 

FOR 

0.  8 

-  1 .  1  M  ICROMETER 

3 

FOR 

3.  5 

-4.0  MICROMETER 

4 

FOR 

8.  5 

-  12.  0  MICROMETER 

5 

FOR 

2  100 

-  3200  MICROMETER 

SRCCOR  A  SINGLY  DIMENSIONED  ARRAY  CONTAINING  THE  X  AND 

Y  COORDINATES  OF  THE  DETONATION  POSITION. 

VALID  RANGE:  -9999.9  -  9999.99  M. 

LOSTRN  A  LOGICAL  VARIABLE  WHICH  IS  .TRUE.  IF  THE 

TRANSMITTANCE  ALONC  A  L I N E-0 F-S I G H T  IS  DESIRED 

AND  .FALSE.  IF  NOT.  IF  .TRUE.  THEN  TRNCOR  AND  RECCOR 

MUST  BE  SPECIFIED. 


C 


TRNCOR  A  SINGLY  DIMENSIONED  ARRAY  CONTAINING  THE  THREE 

COORDINATES  OF  THE  TRANSMITTER.  THE  COORDINATE  SYSTEM 
MUST  BE  IN  METERS.  THE  THIRD  COORDINATE  IS  RESTRICTED 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


TO  BE  BETWEEN  1  AND  5  METERS  (HEIGHT). 
VALID  RANGE:  -9999.9  -  9999.99  M. 


RECCOR 


A  SINGLY  DIMENSIONED  ARRAY  CONTAINING 
COORDINATES  OF  THE  RECEIVER.  (METERS 
THE  THIRD  COORDINATE  MUST  BE  THE  SAME 


THE  THREE 

) 

AS  THE  THIRD 


COORDINATE  OF  TRNCOR. 

VALID  RANGE:  -9999.9  -  9999.99  M. 


EDGE  A  LOGICAL  VARIABLE  WHICH  .TRUE.  IF  CLOUD  DIMENSIONS 
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ARE  TO  BE  CALCULATED  AND  .FALSE.  IF  NOT 
THEN  OBSCOR  AND  SPCHT  MUST  BE  SPECIFIED 


IF  .TRUE 


C 

t 

C 

C  OBSCOR 

C 
C 
C 

C  SPCHT 

C 
c 
c 

C  NEWTIM 

C 

C 

C 

C 

C  TIME 

C 

c 

c 

c 

c 

c 


A  SINGLY  DIMENSIONED  ARRAY  CONTAINING  THE  X  AND  Y 
COORDINATES,  RESP. ,  OF  THE  OBSERVER.  (METERS) 

VALID  RANGE:  -9999.9  -  9999.99 

A  SPECIFIED  HEIGHT  IN  METERS  AT  WHICH  THE  WIDTH  OF 
THE  CLOUD  AS  VIEWED  FROM  POSITION  OBSCOR  IS  DESIRED. 
MUST  BE  BETWEEN  1  AND  5  METERS. 

A  LOGICAL  VARIABLE  WHICH  IS  .TRUE.  IF  ONE  WISHES  TO 
ADVANCE  THE  TIME  AND  FALSE  IF  MORE  RESULTS  ARE 
REQUIRED  AT  THE  CURRENT  TIME.  IF  .TRUE.  THEN  TIME 
MUST  BE  SPECIFIED. 

TIME  MEASURED  IN  SECONDS  AFTER  DETONATION. 

SUCCESSIVE  CALLS  OF  DRTRAN  MAY  KEEP  THE  TIME  UNCHANGED 
FROM  PREVIOUS  CALL  (NEWTIM  -  .FALSE.)  OR  SPECIFY  A 
NEW  TIME  GREATER  THAN  IN  THE  PREVIOUS  CALL.  WHEN 
A  NEW  SOURCE  IS  SPECIFIED  (N EWSRC-. TRUE. )  A  NEW  TIME 
MUST  ALSO  BE  SPECIFIED  (N EWTIM-. TRUE. )  AND  ONLY  IN 
THIS  CASE  MAY  A  TIME  LESS  THAN  THE  PREVIOUS  CALL 
BE  SPECIFIED. 

VALID  RANGE:  OS  -  1000.0  SECONDS. 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 

C 


OUTPUT 

TRNLOS 

I  ERR 

CNTRD 


THE  TRANSMITTANCE  ALONG  THE  L I N E-0 F-S I GHT  BETWEEN 
THE  TRANSMITTER  AND  THE  RECEIVER. 

INTEGER  ERROR  CODE  WHICH  EQUALS  1  IF  A  FATAL  ERROR 
OCCURS  AND  0  OTHERWISE 

A  SINGLY  DIMENSIONED  ARRAY  CONTAINING  THE  HORIZONTAL 
COORDINATE  AND  THE  VERTICAL  COORDINATE  OF  THE 
CENTROID  OF  THE  CLOUD. 


c 

HEIGHT 

THE 

HEIGHT 

OF 

THE  CLOUD 

IN 

METERS 

• 

c 

CENWTH 

THE 

WIDTH 

OF 

THE  CLOUD 

IN 

METERS 

AT 

THE 

CENTROID 

c 

HEIGHT 

W 

c 

SPCWTH 

THE 

WIDTH 

OF 

THE  CLOUD 

IN 

METERS 

AT 

THE 

SPECIFIED 

HEIGHT 
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1 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NCPTS 


THE  NUMBER  OF  POINTS  DETERMINED  ON  THE  EDGE  OF  THE 
CLOUD. 


CPTS 


A  DOUBLY  DIMENSIONED  ARRAY  CONTAINING  THE  COORDINATES 
OF  POINTS  ON  THE  EDGE  OF  THE  CLOUD.  CPTS(1,J) 


IS  THE  HORIZONTAL  COORDINATE 
AND  CPTS(2,J)IS  THE  VERTICAL 
J-TH  POINT.  THE  FIRST  INDEX 


OF  THE  J-TH  POINT 
COORDINATE  OF  THE 
MUST  BE  DIMENSIONED 


TO  2. 


NERR  AN  INTEGER  ERROR  CODE  WITH  VALUES 

-1  CPTS  NOT  LARGE  ENOUCH  FOR  THE  NUMBER  OF 
POINTS  CALCULATED 
0  NO  ERRORS 

1  NO  CONTOUR  FOUNL  BY  THE  CONTOUR  TRACING 
ROUTINE.  THIS  MEANS  THAT  THE  CLOUD  HAS 
DISSIPATED. 

2  CONVERGENCE  WAS  NOT  ACHIEVED  IN  SEARCHING 
FOR  A  POINT  ON  THE  CONTOUR.  MAY  REFER  TO 
A  DISCONTINUITY  IN  THE  CLOUD  EDGE. 

3  CONTOUR  TOO  SMALL  TO  TRACE  I.E.  LESS  THAN 
A  METER  IN  DIAMETER. 

4  LOSTRN  AND  EDGE  WERE  BOTH  SPECIFIED  FALSE 
SO  NO  RESULTS  WERE  CALCULATED 

5  NEWSRC  WAS  FALSE  WHEN  NEWATM  WAS  TRUE. 

A  NEW  ATMOSPHERE  CANNOT  BE  SPECIFIED  DURING 
THE  EVOLUTION  OF  ONE  DUST  CLOUD. 

6  THE  SPECIFIED  HEIGHT  FOR  THE  WIDTH  OF  THE 
CLOUD  WAS  ABOVE  THE  CLOUD 

7  THE  CALCULATION  OF  THE  ATMOSPHERIC 
PARAMETERS  DID  NOT  CONVERGE.  THE  MEASURED 
DATA  MAY  BE  INCONSISTENT. 


SUBROUTINES  CALLED 

ATM  CAL  ACCEPTS  METEOROLOGICAL  DATA  AS  ARGUMENTS  AND  COMPUTES 

NECESSARY  PARAMETERS  IN  COMMON  /WNDPRM/  AND  /TMPPRM/ 

SOURCE  ACCEPTS  SOIL,  CHARGE,  AND  WAVELENGTH  SPECIFICATIONS 

AS  INPUT  AND  COMPUTES  NECESSARY  PARAMETERS  AND  INITIAL 
VALUES  IN  COMMON  /PRT1NF /  AND  /BUOYCL/ 

RISE  GIVEN  CLOUD  DIMENSIONS  DURING  BUOYANT  RISE  DEVELOPMENT 

OF  CLOUD,  RISE  CALCULATES  THE  DIMENSIONS  AT  A  LATER 
TIME 
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C  CONBYN  CONVERTS  CLOUD  DIMENSIONS  AS  SPECIFIED  BY  VARIABLES 

C  USED  IN  RISE  TO  OUTPUT  VARIABLES  DESCRIBING  CLOUD 

C  DIMENSIONS 

C 

C  WIND  COMPUTES  WIND  VELOCITY  AS  A  FUNCTION  OF  HEIGHT 

C 

C  SCON F  DETERMINES  A  SET  OF  POINTS  DEFINING  THE  OBSERVABLE 

C  CONTOUR  OF  THE  CLOUD  AS  VIEWED  FROM  A  DISTANCE  AT  A 

C  CIVEN  ANCLE  TO  THE  WIND  VELOCITY. 

C 

C  PTSPEC  USES  THE  SET  OF  COORDINATE  POINTS  DETERMINED  BY  SCONF 

C  TO  CALCULATE  OUTPUT  VARIABLES  DESCRIBING  CLOUD 

C  DIMENSIONS 

C 
C 

C  FUNCTIONS  CALLED 

C 

C  CWIND  COMPUTES  THE  OPTICALLY  WEIGHTED  CONCENTRATION 

C  FOR  TIMES  AFTER  THE  TRANSITION  FROM  BUOYANT  RISE 

C  TO  WIND  BLOWN 

********************************************************************* 
IERR-0 

CLMRED-. FALSE. 

WRITE(IOOUT, 900) 

900  FORMAT ( 4  9H 1  DIRTRAN-I  DUST  CLOUD  INFRARED  TRANSMISSION  , 

1  12HCALCULAT10N  // 

2  46H  ALL  UNITS  ARE  MKS  UNLESS  OTHERWISE  SPECIFIED.//) 

1 00  R  EAD ( I  0 IN , 701  )N EWATM  ,NEWSRC,LOSTRN, EDGE, NEWT IM 

IF(. NOT. NEW ATM)  GO  TO  10 
IF(ICLMAT. EQ. 1. AND. .NOT. CLMRED)GO  TO  5 
R EAD ( I  0 IN ,  702  )N ATMOS, (ZTMP(I ) , TMPMES ( I ) , ZWND(I ) , 

+W NDMES ( I ), 1-1, 2),ZINV,  THWND 
GO  TO  6 

1PASCT  PASQUILL  CATEGORY 

WNDVEL  WIND  VELOCITY  IN  M/S  MEASURED  AT  2  METERS  ABOVE  GROUND 
WINDIR  WIND  DIRECTION  IN  DECREES  CLOCKWISE  FROM  TRUE  NORTH 

TEMP  TEMPERATURE  IN  DEGREES  C  MEASURED  AT  2  METERS  HEIGHT 

5  NATMOS-I PASCT 

Z WN D  ( 1  )«2. 

ZTMP (1 )-2. 

WNDMES  (1  )  -WNDVEL 
TMPMES  (1  )  -T EMP  +  2  7  3.0 

NOTE  THAT  POSITIVE  X  AXIS  IS  ASSUMED  TO  BE  EAST. 

NOTE  THAT  WINDIR  IS  THE  DIRECTION  OF  ORIGIN  OF  THE  WIND. 
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THWND-2  70.  -WIND1R 
ZINV-150.0 
CLMRED-. TRUE. 

6  CONTINUE 
950  FORMAT (IX) 

WRITE (IOOUT, 950) 

WRITECIOOUT, 951 )N ATMOS 
NIO-1 

I F(NATMOS. EQ. 0 )N 10-2 

WRITE (IOOUT, 952) (ZTMP (I ) , TMPMES (1 ) , ZWND(I ) , WNDMES (I ) , I-l, NIO) 
WRITE (IOOUT, 95 3 )THWND 

951  FORMAT (3  OH  PASQUILL  CATECORY  ,12) 

952  FORMAT  ( 8H  HT  ,F7.2,7H  TEMP  ,F7.2,7H  HT,F7.2,7H  WIND  , 

1  F7.2) 

953  FORMAT (2  2H  WIND  DIRECTION  ,F7 .2) 

10  CONTINUE 


IF ( . NOT. NEWSRC)GO  TO  20 

READ ( 1 0 IN ,  70 3 )N SOIL, NC HRG , CHWT , D ETDE P , D SOD , SRCC0R(1 ) , 
1  S  RCCOR ( 2 ) 

WRITE(IOOUT, 950) 

WRITE (IOOUT, 955 )N SOIL 

WRITE (IOOUT, 954 )NCHRC 

WRITE (IOOUT, 956 )CHWT 

WRITE (IOOUT, 957 )DETDEP 

WRITE (IOOUT, 958 )D SOD 

WRITE (IOOUT, 9 59) (S RCCOR (I ) , I- 1 ,  2 ) 

WRITE(IOOUT,950) 

954  F ORM AT (  1  5H  CHARGE  INDEX  ,12) 


955  F ORMAT ( 1 5H  SOIL  INDEX  ,12) 

964  FORMAT ( 3 OH  WAVELENGTH 

956  F ORMAT ( 3 OH  WEIGHT  OF  CHARGE 

957  F  ORMAT ( 3  OH  DETONATION  DEPTH 

958  F ORMAT ( 3 OH  DEPTH  OF  SOD 

959  F ORMAT ( 3  OH  SOURCE  COORDINATES 
20  CONTINUE 


, F 7 . 2 ,  13H  MICROMETERS 
,F7.  2) 

»  F  7 .  2  ) 

,F7.  2) 

,  2F  1  0.  2  ) 


) 


IF ( .NOT. LOSTRN)GO  TO  30 

READ  (10 IN, 704) (TRNC0R(I ) , I-l, 3), (RECC0R(J ) , J-l, 3) 


I F ( WAVE  1 . 

LT. 

,0. 

4 

)G  0 

TO 

2 

9 

I F ( WAVE  1 . 

CT. 

0. 

7 

)G  0 

TO 

2 

1 

NWL-  l 

GO  TO  30 

I  F  (WAVE  1. 

LT. 

,0. 

8 

)G  0 

TO 

2 

9 

IF  (WAVE  1. 

.  CT. 

,  1. 

1 

)G  0 

TO 

2 

2 

NWL-  2 

GO  TO  30 

I  F  (WAVE  1. 

,  LT. 

.  3. 

5 

)G  0 

TO 

2 

9 

I F (WAVE  1 . 

GT, 

.4. 

0 

)CO 

TO 

2 

3 

29 


n  n  n  n 


NWL-3 
GO  TO  30 

23  IF(WAVE1.  LT.8.  5)GO  TO  29 
IF(WAVE1.GT.  12.0)GO  TO  24 
NWL-4 

GO  TO  30 

24  IF(WAVE  1.  LT.  2  100)GO  TO  29 
IF(WAVE1.GT.3200)GO  TO  29 
NWL-5 

GO  TO  30 

29  WR1TE(I00UT,25) 

25  FORMAT (4  2H  ***  D IRTRAN-I  ERROR  -  WAVE1  OUT  OF  RANGE  ) 

IERR-1 

CO  TO  999 

30  IF  (EDGE)  READ  (IOIN.704)  (0 BSCOR (1 0 ) , 1 0- 1 , 2 ) , S PC HT 
TIMTST-TIME 

IF  (NEWT1M)  READ  (IOIN.704)  TIME 

IF(. NOT. NEWSRC. AND. TIME. LT.TIMTST)GO  TO  997 

WRITE (IOOUT, 960JTIME 

960  F ORMAT ( / /3 OH  TIME  AFTER  BLAST  ,F7.2) 

NERR-0 

CALL  DUSTCL  (N EWATM ,N ATMOS , ZTMP , TMPM ES , ZWN D, WNDM ES , ZI NV, THWK D, 

1  NEWSRC, CHWT,NCHRC,DETDEP, NSOIL, DSOD, NWL, SRCCOR, 

2  LOSTRN, TRNCOR.R ECCOR, EDGE, OBSCOR, S PC HT, NEWT IM .TIME, 

3  TRNLOS, CNTRD, HEIGHT, CENWTH, S PC WTH, NCPT S , CPTS , N ERR ) 

IF  (NERR.EQ.  0)GO  TO  600 

WRITE(IOOUT, 969JNERR 
969  FORMAT ( 2  OH  *****  ERROR  NUMBER  ,12) 

GO  TO  800 

600  IF( .NOT. LOSTRN)GO  TO  700 
WRITE (IOOUT, 950) 

WRITE (IOOUT, 9 64) WAVE l 

WRITE  (IOOUT,  961  )  (TRNCOR  (I), 1*1,3) 

WRITE (IOOUT, 962 ) (R EC COR(I), 1-1,3) 

961  FORMAT ( 3  OH  TRANSMITTER  COORDINATES  , 3F 1 0.  2 ) 

962  F  ORMAT ( 3  OH  RECEIVER  COORDINATES  ,3F10.2) 

WRITE (IOOUT, 970 )TRNLOS 

970  FORMAT ( 4  OH  TRANSMITTANCE  ALONG  THE  LINE-OF-S IGHT  ,E10.  3) 

700  IF( .NOT. EDGE)GO  TO  800 

WRITE (IOOUT, 950) 

WRITE (IOOUT,  96  3) (OBSCOR  (I ),  I- 1, 2) 

963  FORMAT ( 3  OH  OBSERVER  COORDINATES  .2F10.2) 

WRITE (IOOUT ,971  )H EIGHT 

971  FORMAT (2  8H  THE  HEIGHT  OF  THE  CLOUD  IS  , 1  OX, F 1 0. 2 , 7H  METERS) 

C 

990  FORMAT ( ( 1 X, 2F 1 0. 2 )  ) 

C 

WRITE (IOOUT,  9  72) (CNTRD ( I  0 ) , I  0- 1 , 2 ) 
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972  FORMAT (3 OH  THE  CEHTROID  COORDINATES  ARE  ,8X,2F10.2) 

WRITE (IOOUTt973)CENWTH 

973  FORMAT (30H  THE  WIDTH  AT  THE  CENTROID  IS  ,8XF10.  2,  7H  METERS) 

WRITE  (IOOUT  ,  974)SPCHT,  SPCWTH 

974  FORMAT  (1 4H  THE  WIDTH  AT  ,F8.2,11H  METERS  IS  ,5X,F10.2,7H  METERS) 
WRITE (IOOUT, 975 )NCPTS 

975  FORMAT (1 X, I  3, 3 7H  CONTOUR  POINTS  HAVE  BEEN  DETERMINED  ) 

WRITE (IOOUT,  7  08) ( (CPTS(10, I PT ) , 1 0- 1 , 2 ) , I PT- 1 , NCPTS ) 

800  READ(IOIN, 701 )CNTNU 
IF ( CNTNU )C 0  TO  100 

701  FORMAT  (5L1) 

702  FORMAT  (I  1 , 1 X, 1  OF  7 . 2 ) 

703  FORMAT  (2(11, IX),  5F  7. 2) 

704  FORMAT  (6F7.2) 

705  FORMAT  (12) 

708  FORMAT  (  ( 1  X,  2  ( F  1  0.  3 , 2 X )  )  ) 

709  FORMAT (2(12,1X),2F7.2) 

GO  TO  999 

997  WRITE(I00UT,998) 

998  F ORMAT ( 5 4H  ***  DIRTRAN  ERROR  -  TIMES  ARE  NOT  IN  INCREASING  ORDER) 
IERR-1 

999  RETURN 
END 

CONTROLLING  ROUTINE  FOR  D IRTRAN-I  CODE 

********************  SUBFILE  2  *************************** 

SUBROUTINE  DUSTCL (NEWATM ,N ATMOS , ZTMP , TMPM ES , ZWND , WNDM ES , H T I N V , 

1  THWND, NEWSRC.CHWT, NCHRC, DETDEP, N SOIL, D SOD, NWL , SRCCOR, 

2  LOSTRN, TRNCOR, RECCOR, EDGE , OBSCOR, SPCHT, NEWTIM , 

3  TIME, TRNLOS, CNTRD, HEIGHT , CENWTH, S PC WTH , NCPTS , CPTS , N ERR ) 

IMPLICIT  INTEGER*4  (I-N) 

LOGICAL  NEWATM ,N EWS RC , LOSTRN , EDGE , NEWT IM , HOR IZ , ERR 

DIMENSION  ZTMP (2 ) , TMPMES (2 ) , ZWND ( 2 ) , WNDM ES ( 2 ) , SRCCOR ( 2 ) ,TRNCOR(3) 

1  .RECCOR (3), CPTS (2, 200),ORIG(2 ) , TRNFRM ( 2 , 2 ) , TRN ( 3 ) , R EC ( 3 ) 

2  , PNT(3)  .DELS (3 ) , CNTRD ( 2 ) , OB S COR ( 2 ) ,DIR(2) , SCRN(2 ) 

DIMENSION  OWF ( 5 , 2) 

COMMON  /GEOM/COSTH2, SINTH, S INTH2 , V I SEXT , RTP I 
COMMON  /MODE/  HORIZ 

COMMON  /WNDPRM/  D XZ 0 ,  D YXO, DZ 0, U 0, UM , DN , Z I N V 
COMMON  /CLOCK/  FTIME, TWIND 

COMMON  /SEPRTN /  S EP 1 (2 ) , SEP2 (2 ) , PRSEP 1 , PR SEP 2, NUM 1 , NUM 2 
COMMON  /I OUN IT /  I 0 IN , I OOUT , I S PT PF , LOUN IT , ND I RTU , NB SCAT 
EXTERNAL  FUNCT 

DATA  PI/3.  141592  7/ .VISE XT, ZMIN, RES, TANT/ . 1,0. ,1.  ,.1/ 

DATA  ONEM/-1. /.RTPl/l.  7  7245/ 

DATA  OWF/ 1 .  ,  .93,  .52,  .44,  7.  E  - 4,  1.  ,  .93,  .66,  .52,  7.  E  - 4 / 
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C  DAVID  DVORE 

C  AERODYNE  RESEARCH,  INC. 

C  (617)275-9400  X127 

C 

C  PURPOSE 

C 

C  DUSTCL  CALCULATES  DUST  CLOUD  DIMENSIONS  AND  TRANSMITTANCES 

C  THROUGH  DUST  CLOUDS  FOR  GIVEN  METEOROLOGICAL  DATA,  SOIL  TYPE, 

C  BOMB  CHARACTERISTICS,  AND  WAVELENGTH. 

C 

C  SEE  COMMENTS  IN  DRTRAN  FOR  DETAILS. 

C 

C  ******************************************************************* 

IF(LOSTRN.OR. EDGE )GO  TO  100 

NERR-4 

GO  TO  999 

100  I F ( . NOT . NEW ATM )  GO  TO  200 
IF(NEWSRC)  GO  TO  150 
NERR  -5 
GO  TO  999 
150  CONTINUE 

ZINV-HTINV 

CALL  ATMCAL (N ATMOS , ZTMP , TMPM ES , ZWN D , WNDM ES  ,  ERR ) 

IF ( .NOT. ERR)GO  TO  155 
NERR*  7 
GO  TO  999 
155  CONTINUE 
C 

C  COMPUTE  THE  ROTATION  TRANSFORMATION  MATRIX  TO  CONVERT  USER 

C  DEFINED  COORDINATES  INTO  LOCAL  COORDINATES  WITH  X  AXIS  IN 

C  THE  WIND  DIRECTION. 

C 

THETAX-THWND*PI/1  80. 

TRNFRMO  ,  1  )-COS  (THETAX) 

TRNFRM  (  2  ,  2  ) -T  RN  F  RM  ( 1 ,  1  ) 

TRNFRMO  ,  2  )-S  IN  (THETAX) 

TRNFRMO,  1  ) --TRNFRMO,  2) 

200  CONTINUE 

IF (.NOT. NEW SRC)  GO  TO  300 

TWIND-1. E5 

TPRES-0. 

DEL*. 001 

CALL  SOURCE (CHWT.NCHRG, DETDEP, NSOIL.DSOD) 

DO  250  1-1,2 
ORIG (I ) -S  RCCOR (I ) 

250  CONTINUE 

OVRLAP-SPACNC/2. 

SEP  1(1 )-S PACNC*COS(T HARRY -THETAX) 
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SEP1 (2 )-S PACNG*S IN (THARRY-T HE TAX ) 

CALL  PERP (SEP1, SEP2) 

300  CONTINUE 

IF ( . NOT. LOSTRN )  CO  TO  400 

CONVERT  TRNCOR  AND  RECCOR  TO  LOCAL  COORDINATES  WITH  ORICIN  AT 
DETONATION  POINT  AND  X  AXIS  IN  WIND  DIRECTION. 

TRN(3)-TRNCOR(3) 

REC(3)-RECCOR(3) 

DO  320  1-1,2 
TRN  (I  )  -0. 

REC ( I )«0. 

DO  310  J-1,2 

TRN  (I  ) -TRN  (I  )+TRNFRM(l  ,  J)  *  (TRNCOR  (J  )-ORIC(J  )  ) 

REC(1  )-REC(l  )+TRNFRM(I  ,  J)  *  (RECCOR  (J)-ORIC(J)) 

310  CONTINUE 
320  CONTINUE 
400  CONTINUE 

1F(. NOT. EDGE)  GO  TO  500 

COMPUTE  SIN  AND  COS  OF  THE  ANGLE  BETWEEN  THE  OBSERVERS  VIEWINC 
ANGLE  AND  THE  LOCAL  POSITIVE  X  AXIS  WHICH  IS  IN  THE  WIND 
C  DIRECTION. 

C 

CALL  VSUM(ORIG,OBSCOR,ONEM,DIR) 

CALL  UNIT(DIR, DIR, RANGE) 

COSTH-O. 

SINTH-0. 

DO  410  J-1,2 

COSTH-COSTH+T  RNFRM ( 1 , J ) *D IR (J ) 

S  INTH-S  INTH+TRNFRM(2,  J  )  *D  IR  (J  ) 

410  CONTINUE 

SINTH2-SINTH*SINTH 
COSTH2-COSTH**2 
SCRN (1 )-S INTH 
SCRN (2 ) --COSTH 
500  CONTINUE 

IF  (NEWTIM)  CALL  R  IS  E  (T  PR  ES  ,  T  IM  E,  D  EL  ) 

600  IF( . NOT. EDGE  )  CO  TO  650 
FTIME-T  1ME 

CALL  CLDIM (CNTRD, HEIGHT, CENWTH, SPCHT, S PC WTH , NCPT S , CPT S , 

1  NERR) 

I F ( . NOT . ERR )G 0  TO  650 
NERR-6 
GO  TO  999 
650  CONTINUE 

IF (. NOT. LOSTRN )GO  TO  999 
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c 

C  DETERMINE  THE  OPTICALLY  WEIGHTED  CL  VALUE 

C  ALONG  THE  LINE  CONNECTING  THE  TRANSMITTER  AND  RECI EVER 

C 
C 

CALL  VSUM (REC.TRN, ONEM.DIR  ) 

CALL  UNIT(DIR, DIR, RANGE) 

COSTH-D IR  (1 ) 

S  INTH-D  IR  (2  ) 

S  INTH2-S I  NTH  *S I  NTH 
COSTH2-COSTH**2 
SCRN(1 )-S 1  NT  H 
SCRN  (  2 ) «-C  OS  TH 
X-DOTPRD(SCRN  ,  TRN  ) 

HORIZ-. TRUE. 

TRNLOS«EXP(-CWIND(X, Y, TRN (3 ) , TIME) *OWF(NWL, NSOIL) ) 

999  RETURN 
END 

C  CALCULATION  OF  ATMOSPHERIC  PARAMETERS  FOR  DIRTRAN-I  CODE 

C 

C  *****************  SUBFILE  3  ********************* 

C 

SUB  ROUTINE  ATM C AL (N ATM , ZT , TM ES , ZU , UM ES , ERR ) 

IMPLICIT  INTEGER**,  (I-N) 

REAL  M.N 
LOGICAL  ERR 

DIMENSION  ZT (2 ) . TMES(2 ) , ZU(2 ) , UMES(2 ) ,ZL0(6 ) 

COMMON  /WNDPRM/  D XZ 0 , D YXO, DZ 0 , U 0, M , N , Z I N V 
COMMON  /TMPPRM/TO, T 1, TM 

COMMON  /I OUNI T/  I  0  I N , I OOUT , I S PT PF , LOUN I T , N D I RT U , NB S C AT 
DATA  ZL0/-2.  5,  -4.  5,  -  I  3.  5,  1000.  ,55.  ,20.  / 

C  ************************************************************* 

c: 
c 

C  PURPOSE 

C 

C  TU  FIT  THE  BEST  POWER-LAW  PROFILES  OF  W1NDSPEED, 

C  DIFFUSIV1TY,  AND  TEMPERATURE  CONSISTENT  WITH  KNOWN  RELATIONS 

C  GOVERNING  THE  CONSTANT  SHEAR  STRESS  LAYER  TO  GIVEN 

C  MEASUREMENTS  AT  TWO  HEICHTS. 

C 

C 

C  INPUTS 

C 

C  NATM  INTEGER  WHICH  IS  0  IF  WINDSPEED  AND  TEMPERATURE 

C  ARE  AVAILABLE  AT  TWO  HEIGHTS  AND  EQUAL  TO  THE 

C  PASQUILL  CATEGORY  OTHERWISE. 

C 
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SINGLY  DIMENSIONED  ARRAY  CONTAINING  TWO  HEIGHTS 
(IN  METERS)  AT  WHICH  TEMPERATURES  WILL  BE  GIVEN 
MUST  BE  IN  ASCENDING  ORDER. 


C 
C 
C 
C 

C  TMES 

C 
C 
C 

ZU 


UMES 


ZINV 


OUTPUTS 

ERR 


DXZ  0 


DYXO 


DZO 


UO 


M 


N 


TO,  T  1 


SINGLY  DIMENSIONED  ARRAY  CONTAINING  THE  TWO 
TEMPERATURE  MEASUREMENTS  IN  DEGREES  KELVIN 
AT  HEIGHTS  ZT. 

SINGLY  DIMENSIONED  ARRAY  CONTAINING  ONE  OR  TWO 
HEIGHTS  (METERS)  AT  WHICH  WIND  SPEEDS  WILL  BE 
GIVEN.  MUST  BE  IN  ASCENDING  ORDER. 

SINGLY  DIMENSIONED  ARRAY  CONTAINING  THE  ONE  OR 
TWO  WIND  SPEED  MEASUREMENTS  (M/S)  AT  HEIGHTS  UMES. 

INVERSION  HEIGHT  IN  METERS. 


A  LOGICAL  WHICH  IS  TRUE  IF  AN  ERROR  IS  INCURED 
DURING  THE  CALCULATION. 

THE  RATIO  OF  THE  DIFFUSIVITY  IN  THE  X  DIRECTION 
TO  THE  DIFFUSIVITY  IN  THE  Z  DIRECTION.  RETURNED 
IN  COMMON  /WNDPRM/ . 

THE  RATIO  OF  THE  DIFFUSIVITY  IN  THE  Y  DIRECTION 
TO  THE  DIFFUSIVITY  IN  THE  X  DIRECTION.  RETURNED 
IN  COMMON  /WNDPRM/. 

THE  COEFI C I  ENT  OF  Z  *  *N  IN  THE  VERTICAL  PROFILE 
OF  VERTICAL  DIFFUSIVITY.  RETURNED  IN  COMMON 
/WNDPRM/. 

THE  COEFIC I  ENT  OF  Z * *M  IN  THE  VERTICAL  PROFILE 
OF  HORIZONTAL  WIND  SPEED.  RETURNED  IN  COMMON 
/WNDPRM/ . 

THE  EXPONENT  OF  Z  IN  THE  HORIZONTAL  WIND  SPEED 
PROFILE.  RETURNED  IN  COMMON  /WNDPRM/. 

THE  EXPONENT  OF  Z  IN  THE  VERTICAL  DIFFUSIVITY 
PROFILE.  RETURNED  IN  COMMON  /WNDPRM/. 

TM  CONSTANTS  FOR  THE  TEMPERATURE  PROFILE  SUCH  THAT 
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C  T-TO+T  1*Z**TM.  RETURNED  IN  COMMON  /TMPPRM/ . 

C 
C 
C 

C  CALLED  FROM  DUSTCL 

C 

C  NEEDED  FUNCTIONS  AND  SUBROUTINES 

C 

C  TMPC AL  CALCULATES  SCALED  TEMPERATURE  PROFILES 

C 

C  WN  DC  AL  CALCULATES  SCALED  WIND  SPEED  PROFILES 

C 

C 

ERR-. FALSE. 

C  '* 

DELTH  IS  THE  DIFFERENCE  IN  POTENTIAL  TEMPERATURE  BETWEEN  THE 
TWO  HEIGHTS  WHERE  TEMPERATURE  IS  GIVEN. 

Z0-0. 02 

1F(NATM.EQ.0)G0  TO  100 


C 

C 

C 

C 

C 

C 


C 

C 

C 

C 


ASSIGN  ATMOSPHERIC  PROFILE  ACCORDING  TO  SPECIFIED  PASQl’ILL 
CATEGORY 

ZO  FRICTION  HEIGHT 

ZL  MOMIN-OBUKOV  LENCTH 

USTAR  THE  FRICTION  VELOCITY 
TSTAR  THE  SCALING  TEMPERATURE 

Z  L -Z  L  0 (N  ATM ) 

N  P-1 F 1 X (S 1CN  (  1 .  , ZL)  ) 

U  ST AR-UMES ( 1  )  /WNDCAL(Z0, ZL,ZU(1  )  ) 

TSTAR -TMES  (I  )  *U  STAR*  *2/1  .  568/ZL 
I  F ( N  ATM -  4)200,  300,  2  1  0 
100  CONTINUE 

USE  ITERATIVE  PROCEDURE  TO  CONVERCE  ON  BEST  ATMOSPHERIC  PR  OF  I LF 
TO  MATCH  DATA  AT  TWO  HEIGHTS 

DELTH  «TMES(2)-TMES(1  )+.0098*(ZT(2)-ZT(l  )  ) 

N  P-S ICN ( 1 .  .DELTH ) 

DELU-UMES (2 ) -UMES ( 1 ) 

ZULOC-ALOC(ZU(2 )  /Z  U ( 1  )  ) 

Z TLOG -A  LOG ( ZT ( 2 )  /ZT  ( 1  )  ) 

USTAR- (UMES (2) -UMES (1  )  )/ZULOG 
TSTAR-D  ELTH/ZTLOG 
Z  L-.  6  38*TMES (1 ) *U  STAR* *2 /TSTAR 
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t 

i 


IF (ABS (ZL) . GE. 1 000. )G0  TO  300  j 

DO  110  ITER-1,  100  1 

USTAR-DELU/ (WNDCAL(20, ZL, ZU(2 ) )-WNDCAL(Z0,  ZL,  ZU(  1  ) ) )  j 

T  STAR-D  ELTH/ (TMPCAL(Z0,ZL,ZT(2))-TMPCAL(Z0,ZL,ZT(1)))  j 

ZLP-ZL 

ZL-.  63  8*TMES  (1  )  *USTAR*  *2 /T  STAR 
IF  (ABS  (  (ZL-ZLP) /ZLP)  .  LT.  .  0  1  )G0  TO  120 
110  CONTINUE 
ERR-. TRUE . 

CO  TO  999 
1  20  CONTINUE 

Z0-EXP(  .  4*  (WNDCAL  (ZO,  ZL,  ZU(  1  )  )-UMES  ( 1  )  /U  S  T  AR  )  )  *Z  0 
1F(NP)2  00,  300,  2  10 

200  CONTINUE  . 

s 

UNSTABLE  ATMOSPHERE  \ 

i 

DZO-3.  1  5*.  4  * '  S  T  A  R  /  A  B  S  (ZL)  *  *  ( 1 .  /  3 .  )  i 

DXZO-2.  f 

M-1./7.  I 

N-4./3.  1 

CO  TO  430  | 

210  CONTINUE  * 

STABLE  ATMOSPHERE  ? 

i 

DXZO-6. 

M-1./7.  • 

N  -  1  . 

DZ0-.4*USTAR/(l.+47./ZL) 

GO  TO  430 
300  CONTINUE 

NEUTRAL  ATMOSPHERE 

DZ0-.  4*USTAR 
DXZ0-5. 

N  P-0 
N  «  1  . 

M- 1 . / 7 .  j 

430  CONTINUE  I 

COMMON  CALCULATION  TO  UNSTABLE,  NEUTRAL,  AND  STABLE  ATMOSPHERES  j 

IF (ZINV. LE. 10. )ZINV«200.  ' 

UO-UMES  ( 1  )  /ZU  (  1  )  **M 
DYXO- 1 . 

1F(N  ATM.  EQ.  0  )U  0-  (U  0+UMES  (2  )  /ZU  (2  )  **M)  12. 
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IF (NATM . EQ. 0 )G0  TO  450 
ZT  (2  )  -2.  *ZT  (  1  ) 

DELTH-TSTAR*  (TMPCAL(Z0,  ZL,  ZT(2  )  )-TMPCAL(Z0,  ZL,  ZT(  1  )  )  ) 
TMES (2 )-TMES<l )+DELTH 
4  50  T  1-DEL TH/  (ZT(2  )  **M-ZT(1  )  **M) 

TM-M 

TO-TMES  (1  )  -T  1*ZT  (1  ) 

DZ0-. 4 *U  STAR 
N-l. 

999  RETURN 
END 


*****************  SUBFILE  4  ****************************** 

FUNCTION  WNDCAL (Z0, ZL, Z) 

IMPLICIT  INTECER*4  (I-N) 


PURPOSE 

TO  CALCULATE  THE  WIND  SPEED,  U /U * ,  SCALED  BY  THE  FRICTION 
VELOCITY  FROM  GIVEN  FRICTION  HEIGHT  AND  MONIN-OBUKHOV  LENGTH  AT  A 
SPECIFIED  HEIGHT. 


INPUTS 

ZO  THE  FRICTION  HEIGHT  IN  METERS. 

ZL  THE  MONIN-OBUKHOV  LENGTH  IN  METERS. 

Z  THE  HEIGHT  AT  WHICH  THE  SCALED  VELOCITY  IS  DESIRED 

IN  METERS 


RETURNS  VELOCITY  SCALED  BY  FRICTION  VELOCITY 


CALLEb  BY  ATMCAL 


LOGICAL  LOW 

PS  IM  (S  )--  (ALOG  (  (S  +  l.  )  **2*  (S  *S  +  1  )  /8.  )  +2.  *  (ATAN  (  1  .  )  -A  TAN  (S  )  )  ) 
PSIMS(Z) --7. *Z 
PHIM(Z)- (1. -16. *Z) **(-. 25) 

PHIM  THE  SHEAR  OF  MOMENTUM 

C  PSIM  THE  UNIVERSAL  FUNCTION  FOR  THE  DEVIATION  FROM 
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PS  IMS 


LOGARITHMIC  WIND  VELOCITY  BOUNDARY  LAYER  PROFILE  IN  AN 
UNSTABLE  ATMOSPHERE 

THE  SAME  AS  PSIM  FOR  A  STABLE  ATMOSPHERE 


IF(ABS (ZL).LE. 1. E3)CO  TO  100 
WNDCAL-ALOC (1. +Z/Z0) 

GO  TO  999 
100  CONTINUE 

P-S IGN ( 1 . ,ZL) 

LOW-. TRUE. 

S-Z/ZL 

IF(S  .  LE.  1.  5.  AND. S.GE.-2. )CO  TO  10 
S-AMIN  1  (S  ,  1 .  5  ) 

S-AMAX1(S , -  2 .  ) 

LOW-. FALSE. 

10  CONTINUE 

IF (P)l 20, 130, 130 
120  S-l. /PHIM(S ) 

P  S I -PS IM  (  S  ) 

GO  TO  52 
130  CONTINUE 

PS  1-P  S  IM  S  (S  ) 

52  CONTINUE 

WNDCAL-ALOG(l. +Z/Z0)-PSI 
IF (LOW)GO  TO  999 
IF(P)  5  3,  53,  54 

5  3  WNDCAL-WNDCAL+.  7  5-.  9  5*  (-Z  L /Z  )  *  *  ( 1 .  /  3  .  ) 
GO  TO  999 

54  WNDCAL-WNDCAL-15.+10.*Z/ZL 
999  WNDCAL-WNDCAL/. 4 
RETURN 
END 

*******************  SUBFILE  5  ********** 

FUNCTION  TMPCAL (ZO, ZL, Z) 

IMPLICIT  I  NT  EGE  R*  4  (I-N) 


C 

C 

C 

C 

C 

C 

C 

C 

C 


PURPOSE 

TO  CALCULATE  THE  POTENTIAL  TEMPERATURE  SCALED  BY  THE  SCAlE 
TEMPERATURE,  T*,  FROM  GIVEN  FRICTION  HEIGHT  AND  MO N I N -0 BU KH 0 V 
LENGTH  AT  A  SPECIFIED  HEIGHT. 
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INPUTS 


C 

C 

C 

C 

C 

C 

C 

c 


ZO  THE  FRICTION  HEIGHT  IN  METERS. 

ZL  THE  MONIN-OBUKHOV  LENGTH  IN  METERS. 

Z  THE  HEIGHT  AT  WHICH  THE  SCALED  VELOCITY  IS  DESIRED 

IN  METERS 


RETURNS  SCALED  TEMPERATURE 


CALLED  BY  ATM  CAL 


LOGICAL  LOW 

PHIM(Z)«(1.  -16.  *Z  )**(-.  25) 
PSIH (S )--2. *ALOG ( (S  *S  +  1. )  /2.  ) 
PSIHS (Z)--l 1. *Z 


PHIM  THE  SHEAR  OF  MOMENTUM 

PSIH  THE  UNIVERSAL  FUNCTION  FOR  DEVIATION  FROM  LOGARITHMIC 

POTENTIAL  TEMPERATURE  PROFILE  IN  THE  BOUNDARY  LAYER 
OF  AN  UNSTABLE  ATMOSPHERE 

PHIHS  THE  SAME  AS  PHIH  EXCEPT  FOR  STABLE  ATMOSPHERE 

IF(ABS (ZL) . LE. 1. E3)GO  TO  100 
TMPCAL-ALOG (1. +Z/Z0) 

GO  TO  999 
100  CONTINUE 

P«SIGN(1. ,ZL) 

LOW-. TRUE. 

S-Z/ZL 

IF (S .LE. 1.  5.  AND. S.GE. -2. )GO  TO  10 
S -AM  I N  1  (S  ,  1.  5  ) 

S  -AM  AX  1  (S  ,  -  2 .  ) 

LOW-. FALSE. 

10  CONTINUE 

IF  (P  )  1  20,  1  30,  1  30 
120  S-l. / P H I M ( S  ) 

PS1-PSIH  (S  ) 

GO  TO  52 
130  CONTINUE' 

PSI-PSIHS  (S  ) 

52  CONTINUE 

TMPCAL-ALOG (1 . +Z/Z0)-PS  I 
IF (LOW)GO  TO  999 
IF(P  )  5  3,  53,  54 


AO 


5  3  TMPCAL-TMPCAL+.  7 5-.  9 5* (-Z L/Z ) **  ( 1 . / 3 . ) 

GO  TO  999 

54  TMPCAL-TMPCAL-15.+10.*Z/ZL 
999  RETURN 
END 

C  SOURCE  CALCULATION  FOR  DIRTRAN-1  CODE 

C 

c  *******************  SUBFILE  6  *********************** 

C 

C 

C 

SUB  ROUT INE  SOURCE (W , NCHRG, DD , N SOI L, D SOD ) 

IMPLICIT  I  NT  EGE  R*  4  (I-N) 

C  ******************************************************* 

C 

c 

C  PURPOSE 

C 

C  TO  CALCULATE  EXPLOSIVE  DUST  SOURCE  TERM  FOR  THE 

C  DIRTRAN-I  CODE 

C 

C  INPUT 

C 

C  W  THE  WEIGHT  OF  THE  CHARGE  IN  LBS.  TNT 

C  DD  DETONATION  DEPTH  IN  METERS 

C  NSOIL  INTEGER  SOIL  INDEX 

C  DSOD  DEPTH  OF  SOD  IN  METERS 

C 

C  OUTPUT  (RETURNED  IN  COMMON  /PRT1NF  /  AND  /Bl'OYCL/  ) 

C 

C  RO 

C  VGRAV 

C 
C 

C  NPRTS 

C 

C  RSPH 

C  DELT 

C 

C  VXSPH 

C  VZSPH 

C  XCMSPH 

C  ZCMSPH 

C  RISTIM 

C 
C 

C  CALLED  BY  DUSTCL 

C 


INITIAL  CLOUD  RADIUS  IN  METERS 

SINGLY  DIMENSIONED  ARRAY  CONTAINING  OPTICALLY  WEICHTED 
AVERAGE  SETTLING  VELOCITIES  FOR  EACH  SIZE  RANCE  IN 
THE  PARTICLE  DISTRIBUTION  (METERS/SEC) 

THE  NUMBER  OF  SIZE  RANCF.S  IN  THE  PARTITIONING  OF  THE 
PARTICLE  SIZE  SPECTRUM 

THE  INITIAL  RADIUS  OF  THE  CLOUD  IN  METERS 
THE  INITIAL  DIFFERENCE  IN  TEMPERATURE  BETWEEN  THE  CLOUD 
AND  SURROUNDINGS  (DEGREES  KELVIN) 

THE  INITIAL  HORIZONTAL  VELOCITY  OF  THE  CLOUD  (M/S) 

THE  INITIAL  VERTICAL  VELOCITY  OF  THE  CLOUD  (M/S) 

INITIAL  HORIZONTAL  POSITION  OF  THE  CLOUD  (METERS) 
INITIAL  HEIGHT  OF  THE  CLOUD  (METERS) 

TIME  LAPSED  SINCE  DETONATION  IN  SECONDS 
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c  ******************************************************************* 

LOGICAL  HOR1Z 

DIMENSION  CR(5,4),CD(5, 4 ) , OWML (3, 4 ) , OUSV(3, 4),PRTTN(4) 

DIMENSION  S (3) ,BURHTR(5 ),WTRAT(5 ) 

C OMMON /PRTINF/  RO , VGRA V ( 3 ) , N PRTS 

COMMON/BUOYCL/  RSPH.DELT, VXS PH, VZ SPH, XCMSPH, ZCMSPH, 

*  SPHNS (3) .RISTIM 
COMMON  /TMPPRM/  T0.T1.TM 

COMMON  /WNDPRM/  D XZ 0 , D YXO, DZ 0 , U 0, UM , DN , Z IN V 
COMMON  /BURST/  ACCEL, TBURST 

COMMON  /GEOM/COSTH2, SINTH, S INTH2 , V I S EXT , RTP I 
COMMON  /MODE/  HOR1Z 

COMMON/D ISCS/NDSCS,TDSC(1 0),XDSC(10),ZDSC(10),R2DSC(I0), 

1  QDSC(10,3) 

C 

C  CR  IS  THE  CRATER  RADIUS  INDEXED  BY  COEFFICIENT  AND  SOIL  TYPE 

DATA  CR/.483,-.  328,  .0319,  .0536,0.  ,.386, -.291,  .054  3,  .057,  1  1*0.  / 

C  CD  IS  THE  CRATER  DIAMETER  INDEXED  BY  COEFFICIENT  AND  SOIL  TYPE 

DATA  CD/. 261,  -.  453,  .0681,  .224,  .055,  .  198,-.  355,  .0582,  .215,  .0586, 

1  10*0./ 

C  OWML  IS  THE  OPTICALLY  WEIGHTED  MASS  LOADING  COEFFICIENT  INDEXED  BY 
C  BIN  SIZE  AND  SOIL  TYPE 

DATA  OWML/1  .  3E4,  2*0.  ,  2  .  6  1  E  4,  8  *0  .  / 

OWSV  IS  THE  OPTICALLY  WEIGHTED  PARTICLE  SETTLINC  VELOCITY  (CM/SEC) 
INDEXED  BY  BIN  SIZE  AND  SOIL  TYPE 
DATA  OWSV/1  2*0.  / 

C  PRTTN  IS  THE  PARTITIONING  RATIO  INDEXED  ON  SOIL  TYPE 
DATA  PRTTN /4  *. 9/ 

C  BURHTR  IS  THE  RATIO  OF  BURST  HEIGHT  TO  INITIAL  RADIUS  AND  WTRAT 

C  IS  THE  FRACTION  OF  THE  TOTAL  WEIGHT  WHICH  IS  EFFECTIVE  IN  THE  CLOUT 

DATA  BURHTR/0.  ,4.  ,2.  ,4.  ,3. /.WTRAT/. 6,1.  ,.8,1.  ,.6/ 

C 

C 

R1STM-0. 

XCMSPH-0. 

NPRTS-  1 

W  3  -  (W*WTRAT (NCHRG ))**. 3333333 

RO-1  .  5  3 5 *W 3 

TAMB-TO+T  1  *R  0  *  *T  M 

DELT-. 5  7  *T  AMB 

RSPH-RO 

Z  CM  SPH -R  0 

VXS  PH-U  0*ZCMS  PH**UM 
BURHT-BURHTR(NCHRC)*RO 
BURVZ-1.  3  *S  QRT  (R  0  ) 

TBURST-. 1 5*R 0 

VZ  S  PH-  2. *BURHT/TBURST-BURVZ 
ACCEL- (BUR VZ-VZ SPH ) /TBURST 
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CLAM-DD/W3 


CALCULATE  CRATER  RADIUS  AND  DEPTH 

RC-CR(l.NSOIL) 

DC-CD(l.NSOIL) 

IF  (CLAM. LT. 1. E-30)  CO  TO  98 
TERM-1. 

DO  100  1-2, 5 
T  ERM-T  ERM*  CLAM 
RC-RC  +  CR(I ,NSOIL)*TERM 
DC-DC  +  CD(1  ,NSOIL)*TERM 
00  CONTINUE 
98  CONTINUE 
RC-RC*W3 

DC-DC* ( W* WTRAT (NCHRC) )**.  3 
GET  CRATER  VOLUME 

VC-(3.  1  415926/3.  )  *  (RC/DC)**2  *  (DC  -  DSOD)**3 

CALCULATE  OPTICALLY  WEIGHTED  PARAMETERS 

NDSCS-MAX0(5,  I F I  X ( 5 . * W 3 /2 . 4 )  ) 

DO  101  L-l.NPRTS 
S(L)-OWML(L,NSOIL)  *  VC 
VGRAV(L )-CWSV(L , NSOIL ) 

SPHNS  (L  )-PRTTN  (N  SOIL  )  *  S  (L  ) 

QDSC(.l.L)  -  (1 . -PRTTN  (N  SOIL  )  )  *  S  (L  )  /F  LOAT  (N  DSC  S  ) 
101  CONTINUE 

DELH-2. *RO/FLOAT (NDSCS  ) 

Z--DELH/2. 

DO  200  1-1, NDSCS 

Z-Z+DELH 

ZDS  C ( I  ) -Z 

DO  201  J-l.NPRTS 

QDSC (I  ,  J  )  -QDSC( 1 , J ) 

201  CONTINUE 

CON-ALOG (QDSC(1 , 1 ) /V  I SEXT/D ELH/ (2.  *R0) /3.  1 4159) 
IF ( CON  .  GT .  1.  )CO  TO  2  10 
D  -  1 . 

GO  TO  230 
210  D-CON 

DO  220  IT-1,  5 

D-  (CON-1.  +ALOC  (D  )  )  *D  /  (D  -  1 ) 

220  CONTINUE 

230  R  2DS C ( 1 ) -4. *R0*R0/D 

TDSC ( I )--DELH*DELH/D / (DZ0*Z* *DN )  / 4. 
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C 


SICZ-DELH*DELH/D 
WRITE(1  ,  7  7  7  )S IGZ 
XDSC ( I  )  “U  0*Z*  *UM  *  TDSC(I) 

C  WRITE (1 , 7 77 )T DSC (I ) , XDSC (I ) ,ZDSC(I ) , R2DSC ( I ) , QDSC ( I , 1 ) 

7  77  FORMAT  (IX,  5  (1  PE  1  0.  3)  ) 

200  CONTINUE 
RETURN 
END 

C  WIND  DISPERSION  OF  DUST  CLOUD  FOR  DIRTRAN-I  CODE 

C 

C  *****************  SUBFILE  7  ****************************** 

C 

FUNCTION  FUNCT(X.Z) 

IMPLICIT  INTEGER*4  (I-N) 

COMMON  /CLOCK/  TIME, TWI ND 

Q  *********************************************** 

c 

c 

C  PURPOSE 

c 

TO  SUPPLY  A  TRANSMITTANCE  FUNCTION  FOR  THE  CONTOUR  TRACING 
ROUTINE  IN  ORDER  TO  DETERMINE  THE  CLOUD  EDCE. 

INPUT 

X  THE  HORIZONTAL  COORDINATE  IN  METERS 
Z  THE  VERTICAL  COORDINATE  IN  METERS 

OUTPUT 

RETURNS  THE  LOG  OF  THE  OPTICALLY  WEIGHTED  CL  VALUE  FOR  THE 
LINE-OF-S IGHT  SPECIFIED  BY  X, Z 

FUNCTIONS  CALLED 

CWIND 

CALLED  BY  CFUN,  CLIMB,  GRAD  2 

C  ********************************************************** 

Y-0. 

C  WR ITE (1 ,  1  )X, Z, TIME 

1  FORMAT  (I  X,  3  (I  PE  10.  2  )  ) 

IF(Z. LE. 0.  )G 0  TO  100 
EXT-CWIND1X, Y, Z.TIME) 

IF<EXT.LE.  1.  E-30)GO  TO  100 
FUNCT-ALOC  (EXT) 

GO  TO  999 
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100  FUNCT--30. 

999  CONTINUE 
C  WRITEO,  2)FUNCT 

2  F0RMAT(1H+, 30X, 1  PE  10. 2) 
RETURN 
END 


*********************  SUBFILE  8  **************************** 

FUNCTION  CWIND(X, Y, Z,T) 

IMPLICIT  INTECER*4  (I-N) 

REAL  M  , N 
LOGICAL  HORIZ 
DIMENSION  REF(2 ) ,REF0(2  ) 

COMMON /PR T INF/  R0, VGRAV(3 ) , NPRTS 
COMMON  /GEOM/COSTH2, SI NTH, S INT H2 , V I S EXT , RTPI 
COMMON  /MODE/  HORIZ 

COMMON  /WNDPRM/DXZO, DYXO, DZ0,U0,M,N,ZINV 

COMMON /D ISCS/NDSCS.TDSC ( 1 0), XD  S  C ( 1 0) , ZDSC(  1  0) ,R2DSC(  1  0), 

1  QDSC( 1 0, 3 ) 

COMMON  /S  EPRTN /  S E P 1 ( 2 ) , SEP  2 ( 2 ) , PRS E P  1 , PR S EP 2 , NUM  1 , NUM 2 

C  ************************************************************ 

C 

c 

C  PURPOSE 

C 

C  TO  COMPUTE  THE  CONCENTRATION  AT  A  POINT  OR  INTEGRATED  ALONG 

C  A  HORIZONTAL  LINE 

C 

C  INPUT 

C 

C  X , Y , Z  COORDINATES  IN  METERS.  IF  LINE  INTEGRAL  IS  DESIRED, 

C  Y  IS  IGNORED  AND  LINE  IS  SPECIFIED  BY  X  AND  Z. 

C 

C  T  THE  TIME  IN  SECONDS  AFTER  DETONATION 

C 

C  OUTPUT 

C 

C  RETURNS  THE  CONCENTRATION  AT  X,Y,Z,T  IF  HORIZ  IS  .FALSE.  AND 

C  THE  LINE  INTEGRAL  OF  CONCENTRATION  IF  HORIZ  IS  .TRUE. 

C 

C  SUBROUTINES  CALLED 

C 

C  MOMENT  COMPUTES  ZERO  ORDER  MOMENT  AND  INTERPOLATES  FROM 

C  TABLE  OF  HIGHER  ORDER  MOMENTS. 

C 

C  CALLED  BY  FUNCT.DUSTCL 
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COMMON  /PRT1NF/ 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


*  *  * 


************************************************ 


RO  INITIAL  RADIUS  OF  THE  CLOUD  IN  METERS 

VGRAV  SINGLY  DIMENSIONED  ARRAY.  VGRAV(I)  IS  THE  OPTICALLY 

WEIGHTED  AVERAGE  SETTLING  VELOCITY  FOR  PARTICLES  IN  THE 
I  SIZE  RANGE 

NPRTS  THE  NUMBER  OF  PARTICLE  SIZE  RANGES 
COMMON  /DISCS/ 


NDSCS  THE  NUMBER  OF  DISC  SOURCES 

TDSC  SINGLY  DIMENSIONED  ARRAY  CONTAINING  THE  TIME  OF  RELEASE 
OF  THE  DISC  SOURCES 

XDSC  SINCLY  DIMENSIONED  ARRAY  CONTAINING  THE  X  COORDINATE 
OF  THE  CENTER  OF  THE  DISC  SOURCES 
ZDSC  SINGLY  DIMENSIONED  ARRAY  CONTAINING  THE  Z  COORDINATE 
OF  THE  CENTER  OF  THE  DISC  SOURCES 
R2DSC  SINGLY  DIMENSIONED  ARRAY  CONTAINING  THE  SQUARE  OF  THF 
RADII  OF  THE  DISC  SOURCES 

QDSC  DOUBLY  DIMENSIONED  ARRAY.  QDSC(I,J)  IS  THE  NUMBER  OF 
PARTICLES  OF  THE  J  SIZE  RANCE  IN  THE  I  DISC. 


SUM  THE  CONTRIBUTIONS  OF  THE  DISC  SOURCES  TO  THE 
OPTICALLY  WEIGHTED  CONCENTRATION  AT  X,Y,Z,T 

CWI ND-O. 

DO  2  11  I  -  1 , NDSCS 
TOF-T-TDSC  ( I  ) 

REF0(1  ) "XDSC ( I  ) 

REFO (2  )  -0. 

ROH  2-R  2DSC  (I  ) 

H-ZDSC  (I  ) 

IF (HOR1Z)  REFO  (  1  ) -R EFO  (  1  ) *S INTH 
DO  210  J-l, NPRTS 

DETERMINE  MOMENTS  FOR  CURRENT  SOURCE  DISC  AT  Z 

CALL  MOMENT (VGRAV(J ), Z.H.TOF ,Q,XBAR, SICW2, SICP2)  f 

I F (Q . GT. 1. E-20)  GO  TO  113 

IF  Q  IS  TOO  SMALL,  ITS  CONTRIBUTION  IS  INGORED 


i 


CWNDSC-0. 
GO  TO  210 
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1  1  3  CONTINUE 

RX2*R  OH  2+2.  *SICW2 
RY2-R0H2+2.  *SICP2 
IF(HORIZ)GO  TO  120 

*****  COMPUTE  CONCENTRATION  AT  (X,Y,Z,T) 

ARC“-(X-REF0(1  )  -XB  AR  )  *  *2  /  R  X  2 
I F  ( AB  S  '  RC).CT.  30.  )CO  TO  150 
CWNDSC  (Q/RTP1/SQRT(RX2))*EXP(ARG) 

ARG*-(Y-R£F0(2)  )**2/RY2 
IF(ABS(ARC) . CT. 30. )GO  TO  150 
C  Y*E  XP (ARC  )/RTPl/SQRT(RY2) 

CWNDS  C“QDSC ( I  , J ) *CWNDS C* C Y 
GO  TO  150 

COMPUTE  CONCENTRATION  ALONC  L I N E -0 F-S 1 C H T  SPECIFIED  BY  X,Z 

120  CONTINUE 

REF  (  1  )  -REFO  (1  ) 

REFF2-RX2*S1NTH2+RY2*C0STH2 

ARC*- (X-REF ( 1 )-XBAR*S INTH ) *  *2 /REFF2 

I  F  (ABS  (ARG  )  .  GT.  30.  )CO  TO  150 

CWNDSC-EXP(ARC)  ‘  /  S  QR  T  ( R  E  F  F  2  )  /R  T  P  I 

CWNDSC-CVNDSC*Q  *QDSC ( 1  ,  J  ) 

150  CONTINUE 

CUIND-CWIND+CWNDSC 
210  CONTINUE 
2 1  1  CONTINUE 
RETURN 
END 

*********************  SUBFILE  9  **************************** 

SUBROUTINE  CONVRT  (T  ) 

IMPLICIT  INTEGER*^  (I-N) 

REAL  M  ,  N 
LOGICAL  HORIZ 

COM  MON/BUOYCL/RSPH,DELT,VXSi  1 , VZ  S PH , XC M S  Ph ,  ZCMSPH, 

1  SPHNS (3  )  .RISTIM 
COMMON/PRTINF/  R0,VGRAV(' ) , KPRTS 
COMMON  /WNDPRM/DXZO, DYXO, DZO, UC , M  ,  N  ,  Z  IN  V 
COMMON  /G  EOM /COSTH2,  SINTH, SINTH2.V1SEXT, RTPI 
COMMON  /MODE/  HORIZ 

COMMON/DISCS/NDSCS,TDSC( 1 0),XDSC( 1 0) , ZDSC( 1 0)  ,R2DSC( 1 0), 

1  QDSC ( 1 0, 3) 

COMMON  /CLOCK/  TIME.TWIND 


oonnooonnoonnooononorioonoooo 


PURPOSE 


C 
C 

c 
c 

C  TO  CONVERT  THE  CURRENT  BUOYANT  CLOUD  INTO  DISC  SOURCES  FOR 

C  USE  BY  THE  WIND  DISPERSION  MODEL 

C 

C  INPUT 

T  THE  TIME  IN  SECONDS  AFTER  DETONATION 

OUTPUT 


SUBROUTINES  CALLED 

MOMENT  COMPUTES  ZERO  ORDER  MOMENT  AND  INTERPOLATES  FROM 
TABLE  OF  HICHER  ORDER  MOMENTS. 

CALLED  BY  FUNCT.DUSTCL 


COMMON  /BUOYCL/ 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


RS  PH  RAD  U I S  OF  THE  SPHERICAL  BUOYANT  CLOUD  IN  METERS 
DELT  TEMPERATURE  EXCESS  OF  CLOUD  OVER  AMBIENT  ATMOSPHERE 
AT  SAME  HEIGHT  IN  DEGREES  KELVIN 
VXSPH  X  COMPONENT  OF  CLOUD  VELOCITY  IN  METERS/SEC 
VZSPH  Z  COMPONENT  OF  CLOUD  VELOCITY  IN  METERS/SEC 
XCMSPH  X  COORDINATE  OF  CENTER  OF  SPHERE  IN  METERS 
ZCMSPH  Z  COORDINATE  OF  CENTER  OF  SPHERE  IN  METERS 
SPHNS  SINGLY  DIMENSIONED  ARRAY.  SPHNS(I)  IS  THE  NUMBER 
OF  PARTICLES  OF  THE  1  SIZE  RANGE  IN  THE  SPHERE. 

RISTIM  THE  TIME  IN  SECONDS  AFTER  DETONATION  THAT  THE  CLOUD 
HAS  RISEN  BUOYANTLY 

COMMON  /PRTINF / 

RO  INITIAL  RADIUS  OF  THE  CLOUD  IN  METERS 

VCRAV  SINGLY  DIMENSIONED  ARRAY.  VCRAV(I)  IS  THE  OPTICALLY 

WEIGHTED  AVERACE  SETTLING  VELOCITY  FOR  PARTICLES  IN  THE 
I  SIZE  RANGE 

NPRTS  THE  NUMBER  OF  PARTICLE  SIZE  RANCES 
COMMON  /DISCS/ 
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c 

NDSCS 

THE  NUMBER  OF  DISC  SOURCES 

c 

c 

TDSC 

SINGLY  DIMENSIONED  ARRAY  CONTAINING 
OF  THE  DISC  SOURCES 

THE 

TIME  OF  RELEASE 

c 

c 

XDSC 

SINGLY  DIMENSIONED  ARRAY  CONTAINING 
OF  THE  CENTER  OF  THE  DISC  SOURCES 

THE 

X  COORDINATE 

c 

c 

ZDSC 

SINGLY  DIMENSIONED  ARRAY  CONTAINING 
OF  THE  CENTER  OF  THE  DISC  SOURCES 

THE 

Z  COORDINATE 

c 

c 

R2DSC 

SINGLY  DIMENSIONED  ARRAY  CONTAINING 
RADII  OF  THE  DISC  SOURCES 

THE 

SQUARE  OF  THE 

c 

QDSC 

DOUBLY  DIMENSIONED  ARRAY.  QDSC(I, 

J)  IS  THE  NUMBER  OF 

c 

c 

NSPH-5 
D  EL*  2 . 

PARTICLES  OF  THE  J  SIZE  RANGE  IN  THE  I 

/FLOAT (NSPH ) 

DISC. 

A2-3.  /FLOAT  (NSPH  )  **2 
A3-2.  /  FLOAT  (NSPH)  *  *  3 
DELZ  *DEL*RSPH 
HREF-ZCMSPH-R  SPH-. 5*D  ELZ 
R  2  *R  S  PH  *R  S  PH 
DO  50  1*1, NSPH 
IDSC-NDSCS+I 

ZD5C(IDSC)*HREF+FLOAT  (1  )*DELZ 
HMZ-ZDSC  (  1DSC)  -Z  CMS  PH 
RD2-R  2-HMZ*HMZ 

VFRAOA2*  (  FLOAT  (2  *1-1  )  )  -A  3  *F  LOAT  (  3  *1  *  ( I  -  1  )  +  1  ) 

DO  5  1  PR  T* 1 , NPRTS 

QDSC(IDSC, IPRT)*VF  RAC  *S  PH  NS (I  PR  T ) 

5  CONTINUE 

A-ALOG  (QDSC(lDSC.l) /V  I S EXT /D ELZ / S QRT (R D 2 )  /RTPI**2 ) 

IF ( A . GT .  1.  ) G  0  TO  210 
D-  1. 

CO  TO  230 
2  10  D*A 

DO  220  IT*  1 , 5 

D-  (A -I  .  +ALOG  (D  )  )  *D  /  (D  -1  ) 

220  CONTINUE 

230  R2DSC(IDSC)«RD2/D 

GAPTIM«-DELZ*DELZ/D/4./(DZ0*ZDSC(IDSC)**N) 

XDSC ( IDSC) -U 0*ZDSC ( I  DSC ) **M  *  CAPTIM  +XCMSPH 
TDSC(IDSC)-T+GAPTIM 
50  CONTINUE 

NDSCS-NDSCS+NSPH 

TWIND-T 

RETURN 

END 

BUOYANTLY  RISING  DUST  CLOUD  FOR  DIRTRAN-I  CODE 
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c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  R  I  S  E  (T  PR  E  S  ,  TN  EXT  ,  D  E  L  ) 

IMPLICIT  INTECER*4  (1-N) 

REAL  M.NDIF 
DIMENSION  W  K  ( 1  2,  6) 

C OMMON /BUOYCL/  Y(6),SPHNS(3),RISTIM 

COMMON  /WNDPRM/  D XZ 0 , D YXO , DZO, U 0, M , ND I F , Z 1 N V 

COMMON  /CLOCK/  TIME.TWIND 

COMMON  /I  OUN  IT  /  1  0  I N  ,  I  OOUT  ,  I  S  PT  PF  ,  LOUN  I  T  ,  N  D  1  R  T  l' ,  B  SC  A  T 
DATA  HMIN.ACCURC.UK, N,ND/.  001,.  001,  72*0.  ,12,  12/ 

******************************************************************* 


PUR  POSE 

THIS  ROUTINE  CALLS  A  RUNCE  KUTTA  ROUTINE  TO  INTEGRATE  IN  TIME 
THE  EQUATIONS  FOR  THE  RISE  OF  A  BUOYANT  CLOUD  BEC1NNINC  AT  TPPE5 
ENDING  AT  TNEXT  UNLESS  THE  CONDITION  FOR  SWITCHING  TO  THE  WIND 
DISPERSION  MODEL  IS  ENCOUNTERED  IN  WHICH  CONVRT  IS  CALLED. 

SEE  SUBROUTINE  DIFEQ  FOR  THE  DEFINITIONS  OF  Y(I). 


ARGUMENTS 

TPRES  AS  INPUT  TPSES  IS  THE  INITIAL  TIME  OF  THIS  SEGMENT  OF 

INTEGRATION  AND  IS  RETURNED  WITH  THE  VALUE  OF  THE  LAST 
SUCCESSFUL  INTEGRATION  STEF. 

TNEXT  THE  ENDTPOINT  OF  THE  TIME  INTERVAL  WHICH  IS  INPUT. 

REQUIRED  SUBROUTINES 

RKM  A  RUKCE-KUTTA-MERSON  INTEGRATION  ROUTINE 

CONVRT  A  SUBROUTINE  WHICH  CONVERTS  THE  CURRENT  BUOYANT 
DUST  CLOUD  TO  A  NUMBER  OF  DISC  SOURCES  FOR  THE 
WIND  DISPERSION  MODEL.  A  GAP  TIME  DURING  WHICH  THE 
BUOYANT  MODEL  IS  CONTINUED  IS  COMPUTED. 

CALLED  BY  DUSTCL 


********** 


************* 


I F (TNE XT . GT. TW IND )GO  TO  999 
T  2-T  PRES 
C 

C  PERFORM  INTEGRATION  IN  SEGMENTS  OF  TIME 

C 

10  DO  20  NT« 1 , 40 
T  1-T  2 
T  2-1. 2*T  1 
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I  F (T  2. LE. 0.  ) T 2*. 5 

IF(T  2. GT. TNEXT )T  2 -T NEXT 

IF(Y  (1  )+Y(6  ).CE.  ZIKV)GO  TO  200 

IF  (DEL.LT.HMIN)DEL-HMIN 

CALL  RKM(N,T1,T2,Y,HM1N,DEL,  ACCURC , WK , ND ) 

C 

C  CHECK  TO  SEE  IF  CLOUD  GROWTH  IS  DOMINATED  BY  WIND  DIFFUSION 

C  OVER  BUOYANT  RISE  BY  COMPARING  WIND  DIFFUSIVITY,  DIFW,  TO 

C  THE  EFFECTIVE  BUOYANT  DIFFUSIVITY,  DIFB. 

C 

5  DIFB-ABS  (.  25*Y  (1  )  *Y  <4  )  ) 

D  IFW-DZO*  (Y  (6  )+Y  ( 1  )  )  **N  DIF 
IF (DIFB. GT. DIF W)CO  TO  15 
CALL  CONVR  T  (T  2  ) 

GO  TO  200 


C 

C 

C 


15  CONTINUE 

IF  (T 2. CE. TNEXT )CO  TO  200 
I  F(T  2.  GT.  300.  )GO  TO  9  9 
20  CONTINUE 
99  WRITE(IOOUT,  98) 

98  FORMAT ( 5  5H  ***  DIRTRAN-I  ERROR  -  5  MINUTE  CUT-OFF  ON  BUOYANT  RISE 
1) 

STOP 

200  TPRES-T  2 

R  1STIM-TPRES 
9  99  RETURN- 
END 


*****  SUBFILE  1  1 


********* 


SUBROUTINE  D  I  F  EQ  (N  ,  T  ,  Y  ,  Y  P  ) 

IMPLICIT  INTEGERS  (I-N) 

DIMENSION  Y  (N  )  ,  YP  (N  ) 

COMMON /PRT INF /ROC  L, VCRAV(3 ) , NPRTS 
COMMON  /TMPPRM/TO,  T  1,  TM 

COM  MON/WNDPRM/DXZO, D YXO , DZ 0 , U 0, UM , DN , 2  I  N V 
COMMON  /BURST/  ACCEL,  TBl’RST 
DATA  ALPHAK/. 2  5/ 

C  ****************************************************************** 

c 

c 

C  PURPOSE 

C 

C  D IF  EQ  CONTAINS  THE  PARTIAL  DIFFERENTIAL  EQUATIONS  FOR  THE 

C  RISE  OF  A  BUOYANT  CLOUD  WHICH  ARE  USED  BY  SUBROUTINE  R KM . 

C 

c 

C  INPUT 
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C 

C  K 

C  T 

C  Y  (  1  ) 
C  Y(2) 

Y  ( 3  ) 

Y  ( 4  ) 

Y  (  5  ) 

Y  ( 6  ) 


THE  NUMBER  OF  DEPENDENT  VARIABLES 
THE  INDEPENDENT  VARIABLE,  I.E.  TIME 
RADIUS  OF  CLOUD 

CLOUD  TEMPERATURE  MINUS  SURROUNDING  TEMPERATURE 
CLOUD  VELOCITY  IN  WIND  DIRECTION- 
VERTICAL  VELOCITY  OF  CLOUD 

X-COORDINATE  OF  CENTER  OF  MASS  FOR  THE  CLOUD 
THE  HEIGHT  OF  THE  CLOUD  C.O.M. 


OUTPUT 

YP  AN  ARRAY  CONTAINING  COMPUTED  DERIVATIVES  OF  THE  DEPENDENT 
VARIABLES  WITH  RESPECT  TO  THE  INDEPENDENT  VARIABLE. 


REQUIRED  FUNCTIONS 
NONE 


CALLED  BY  R KM 


;  *  *  ★  *  A  A 


IF(T.LT.TBURST)CO  TO  200 
DELT-T  1  * Y (6  )  **TM 
TA-T  0+DELT 
DTADZ «TM*DELT/Y  (6  ) 

VW-U  0*Y  (6  )  **UM 

VR-S QRT ( (Y (3 )-VW) *(Y (3 )-VW)+Y(4 ) *Y (4  )  ) 


TA  THE  AMBIENT  ATMOSPHERIC  TEMPERATURE  AT  CLOUD  HEIGHT 

DTADZ  THE  TEMPERATURE  GRADIENT  AT  CLOUD  HEIGHT 

VW  THE  WIND  VELOCITY  AT  CLOUD  HEIGHT 

VR  THE  RELATIVE  VELOCITY  OF  THE  CLOUD  WITH  RESPECT  TO  W! 

TR  THE  RATIO  OF  CLOUD  TEMPERATURE  TO  AMBIENT  TEMPERATURE 


TR-Y (2 )  /TA 


CALCULATE  ARVOL,  THE  SURFACE  AREA  TO  VOLUME  RATIO 
WHICH  DEPENDS  ON  THE  NUMBER  AND  PLACEMENT  OF  CHARCES 

ARVOL-3. / Y ( 1 ) 

10  CONTINUE 

YP  ( 1  )  mAlPHAK*VR 

YP(2  >-  -  (1  .  +T  R)  *ARVOL*Y  (2  )  *Y  P(1  )  -Y  (4  )  *  (DTADZ  ) 
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YP(3)*ARVOL*(VW-Y(3)  )  *Y  P  (  1  ) 

YP(4  )-9.  8*TR-ARV0L*Y  (4  )  *Y  P  (  1  ) 

YP(5)«Y(3) 

YP(6  )-Y(4  ) 

GO  TO  999 
2  00  N  *=  6 

DO  2  10  1*=  1 ,  N 
YP (I )-0. 

210  CONTINUE 

Y  P  (  4  )  -A  C  C  E  L 
YP(5  )=Y  (3  ) 

Y  P  (  6  )  =Y  (  4  ) 

9  99  RETURN 

END 

C 

C  ********************  SUBFILE  12  ****************************** 

C 

SUBROUTINE  CLDIM (CNTRD, HEIGHT, CENWTH, S PC HT,S PC WTH.N OPTS, CRTS  5, 

1  NERR) 

IMPLICIT  I  N  T  EG  E  R*  4  (I-N) 

Q  ******************************************************************* 

c 

c 

C  PURPOSE 

C 

C  CLDIM  CALCULATES  FIVE  CONTOUR  POINTS  AND  CLOUD  DIMENSIONS  AS 

C  SEEN  FROM  THE  SPECIFIED  OBSERVER  POSITION.  CLDIM  REQUIRES  CLOUD 

C  PARAMETERS  FROM  THE  BUOYANT  RISE  STACE  OF  CLO'JD  DEVELOPMENT  WHICH 

C  ARE  SUPPLIED  IN  COMMON  STORAGE  /BUOYCL/  AND  /PRTINF /  AS  WELL  AS 

C  VIEWING  GEOMETRY  WHICH  IS  SUPPLIED  IN  COMMON  /GEOM/.  SPCKT  IS 

C  REQUIRED  INPUT  IN  THE  ARGUMENTS.  ALL  OUTPUTS  ARE  ARGUMENTS. 

C 
C 

C  INPUT 

C 

C  S  PC  1!  T 

C 
C 
C 
C 

C  OUTPUT 

c 

C  CNTRD 

C 
C 

C  HEIGHT 

C  CENWTH 

C  SPCWTH 


THE  SPECIFIED  HEIGHT  AT  WHICH  THE  WIDTH  OF  THE  CLOUD 
IS  DESIRED.  (METERS) 


A  SINGLY  DIMENSIONED  ARRAY  OF  LENGTH  2  WHICH  CONTAINS 
THE  HORIZONTAL  AND  VERTICAL  COORDINATES,  RESP.,  OF  THE 
CLOUD  CENTROID.  (METERS) 

THE  HEIGHT  OF  THE  CLOUD  IN  METERS 

THE  WIDTH  OF  THE  CLOUD  AT  THE  CENTROID  HEIGHT  IN  METERS 
THE  WIDTH  OF  THE  CLOUD  AT  THE  SPECIFIED  HEIGHT  (METERS) 
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NC  PTS 
CPTS 


C 


THE  NUMBER  OF  CONTOUR  POINTS  (-5) 

A  DOUBLY  DIMENSIONED  ARRAY  OF  SIZE  (2 , N ) , N . C E . 5 ,  WHICH 
CONTAINS  THE  HORIZONTAL  AND  VERTICAL  COORDINATES  OF 
THE  FIVE  CONTOUR  POINTS.  (METERS) 


REQUIRED  SUBROUTINES 

SCONF  MAIN  ROUTINE  FOR  CONTOUR  TRACING  ALGORITHM.  USED 
WITH  WIND  MODEL. 

CALLED  BY  DUSTCL 

D IMENS ION  CNTRD (2 ) , C  PTS (2  ,  200),CPTS5(2,5) 

LOGICAL  HORIZ, NOCONT 

COMMON  /BUOYCL/RSPH , DELT, VX,VZ,XCM,ZCM,SPHNS(3),TIM 
COMMON  /PRTINF/RO,  VGRAVO  )  ,NPRTS 
COMMON  /GEOM/COSTH2, SINTH, SINTH2, VISE XT , RTPI 
COMMON  /MODE/  HORIZ 
COMMON  /CLOCK/  T.TW1ND 

COMMON  /WNDPRj./DXZO,  DYXO,  DZO,  U  0,  UM  ,  DN  ,  Z  I  N  V 

COMMON/DISCS/NDSCS,TDSC(10),XDSC(10),ZDSC(10),R2DSC(1C), 

1  QDSC  (  10,3) 

COMMON  /SPECS/  RES , STEP, TANT, CON 
EXTERNAL  FUNCT.GFUN 

DATA  VISEXT.ZM1N, RES, TANT  / .  1 , 0 .  ,  . 4 ,  . 1  / 

HORIZ-. TRUE. 

CON-ALOC (VISEXT  ) 

C  PTS  5(2,  1  ) -S  PCHT 
C  PT  S  5  (  2  ,  5)«SFCHT 
U-U  0*C  PTS  5(2,  1 ) **UM 
CPTS5  ( 1 ,  1  )  -T  *U *S INTH 
CPTS  5  (I  ,  5  )  -C  PTS  5  (1  ,  1  ) 

N  SERCH--I 
S  TEP-  2  0. 

CALL  CLIM B (FUNCT ,GFUN, CPTS 5, FP  I, NSERCH, NOCONT ) 

N  SERCH- 1 
S  T  EP -  2  0. 

CALL  CLIMB(F UN CT, CF UN, CPTS5 ( 1 , 5), FP1, NSERCH, NOCONT) 
SPCWTH«CPTS5(1  ,  5  ) -C  PTS  5(1  ,  1  ) 

NCPTS*=5 

IF (T . LE. TWIND)GO  TO  50 

IND-NDSCS-2 

C  PT  S  (  2 ,  1 ) “Z  DSC ( I ND ) 

C  PTS  ( 1  ,  1)«(XDSC(1ND)  +  (T-TDSC(IND)  )  *U  0*ZDSC  (  I  ND  )  **UM)  *S  INTH 
STEP-20. 

CALL  SCONF ( FUNCT, ZMIN, NCPTS.CPTS.NERR) 


54 


C  I F (N  ERR . N  E. 0)  GO  TO  999 

CALL  PTSPEC(NCPTS,  CPTS,  XH1CH,  HEIGHT,  XC  1  ,  CNTKD ,  XC2,  CENL'TF  ) 

CO  TO  100 

50  CNTKD  (1  ) -XC  M*  S INTH 
CKTRD  (2  )  -Z  CM 
XH1GH*=CNTRD(  1  ) 

H  E 1 G  H  T  =  Z  CM+R S  PH 
C  E  NW  T  H  =  2  .  *  (RSPH  ) 

100  C  PTS5  ( 1  ,  2  )  -CNTRD  ( 1  )  -C  ENWTH/2. 

CPTS5 (1 , 4 )-CNTRD  (  1  )+C  ENWTH/2 . 

CPTS5  ( 1  ,  3  )  =  XH  IC1I 
CPTS5 (2,2) =C  N  TRD  (2  ) 

C  PTS5  (2,  3  )  =H  EICHT 
CPTS5 (2,4) -CNTRD  (2  ) 

N  C  FT  S  =  5 
9  9  9  RETURN 
EN  I) 

C:  ROUTINE  FOR  TRACING  A  CONTOUR  OF  A  FUNCION  OF  TWO  VARIANT 

l 

C  ******************  SUBFILE  13  ********************************* 

( 

< 

C  SCONF  IS  THE  CONTROLLING  ROUTINE  FOR  THE  CONTOUR  TRACING 

C  ALGORITHM.  THE  FUNCTION  WHOSE  CONTOUR  IS  DESIRED  MUST  BE  A 
C  CONTINUOUS,  REAL  VALUED  FUNCTION  OF  TWO  VARIABLES. 

C  THE  CONTOUR  FOUND  IS  ONE  CONTINUOUS  CLOSED  CONTOUR  OR 

C  ONE  CONTINOUS  CURVE  BEGINNING  AND  ENDINC  ON  THE  SPECIFIED 
C  BOUNDARY.  THERE  MAY  BE  OTHER  PIECES  TO  THE  CONTOUR. 

C 

C  INPUT 
C 

C  FUNCT 

C  CON 

C  YM  N 

C  RES 

C  DELTA 

C 

C  THETAN 

C 
C 

C  MAXDIM 

C 

C  CP 

C 
C 
C 
C 


THE  FUNCTION  WHOSE  CONTOUR  IS  TO  BE  FOUND 
THE  CONTOUR  LEVEL. 

LOWER  BOUND  FOR  THE  SECOND  COORDINATE 
THE  RESOLUTION  LENGTH. 

THE  INITIAL  STEP  SIZE  WHEN  LOOKINC  ALONG  THE 
GRADIENT  TO  FIND  A  POINT  ON  THE  CONTOUR. 

THE  TANCENT  OF  THE  MAXIMUM  ALLOWABLE  ROTATION  ANCLE 
BETWEEN  SUCCESSIVE  LINE  SEGMENTS  OF  THE  POLYGONAL 
CONTOUR  . 

THE  NUMBER  OF  POINTS  FOR  WHICH  STORAGE  HAS  BEEN 
ALLOCATED  IN  THE  ARRAY  CP. 

A  DOUBLY  DIMENSIONED  ARRAY  TO  BE  FILLED  WITH  THE 
COORDINATES  OF  THE  CONTOUR  POINTS.  CP(I,J)  IS  THE 
I-TH  COORDINATE,  1-1,2  ,  OF  THE  J -T H  POINT  OF  THE 
CONTOUR.  UPON  CALLING  SCONF,  C P ( I  ,  1  )  , I  -  1  , 2  ,  SHOULD 
BE  THE  COORDINATES  OF  THE  BEST  GUESS  AS  TO  WHERE  THE 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


OUTPUT 

CP 


MAXD IM 
ERROR 


CONTOUR  IS  LOCATED. 


A  DOUBLY  DIMENSIONED  ARRAY  FILLED  WITH  THE 
COORDINATES  OF  THE  CONTOUR  POINTS.  CP(1,J)  IS  THE 
I-TH  COORDINATE,  1-1,2  ,  OF  THE  J -T H  POINT  OF  THE 
CONTOUR . 

THE  NUMBER  OF  POINTS  REPORTED  IN  CP. 

A  NUMERICAL  FLAG  TELLING  WHAT  TYPE  OF  ERROR 
OCCURRED  WHILE  RUNNING,  IF  ANY. 

MEANING  OF  ERROR  CODE- 
-1 =  >A  R  R  A  Y  FILLED  COMPLETELY 
0«>CONTOUR  CLOSED  OR  MET  BOUNDARIES,  NO  PROBLEM. 
1  *  >N  0  CONTOUR  FOUND 

2->PTFIND  UNABLE  TO  FIND  NEXT  POINT.  CHECK  FOR 
KINK  OR  DISCONTINUITY  OF  THE  FUNCTION. 


CALLED  SUBROUTINES. 

CLIMB  LOCATES  FIRST  POINT  OF  CONTOUR  OR  THAT  CONTOUR 
DOESN'T  EXIST. 

PTFIND  FINDS  ADDITIONAL  POINTS  ON  THE  CONTOUR  GIVEN  AT  LEAST 
ONE. 


LOCAL  VARIABLES 


LEFT 

TMP 

NOCONT 

ENDCON 


A  LOGICAL  WHICH  IS  .TRUE.  FOR  CO UNT ER CLOC KW  I  S E  CONTOUR 
SEARCH  AND  .FALSE.  FOR  CLOCKWISE  SEARCH 
TEMPORARY  STORACE  USED  FOR  SWAPPING  POSITIONS  IN 
THE  ARRAY  OF  CONTOUR  POINTS. 

A  LOGICAL  ERROR  FLAG  MEANING  EITHER  THE  CONTOUR 
DOES  NOT  EXIST  OR  IS  TOO  SMALL  TO  BE  TRACED 
A  NUMERICAL  FLAC  INDICATING  THE  STATUS  OF  THE 
CONTOUR  POINTS  TRACED  WITH  VALUES: 

-1  THE  ARRAY  CP  IS  FILLED 

0  THE  CONTOUR  IS  CLOSED 

1  THE  CONTOUR  HAS  MET  A  BOUNDARY 

2  THE  REQUIRED  STEP  SIZE  HAS  BECOME  LESS 
THAN  THE  RESOLUTION 


SUBROUTINE  SCONF ( FUNCT , YM N , MAXD IM , CP, ERROR ) 
IMPLICIT  INTEGER**.  ( I -N  ) 

EXTERNAL  FUNCT, GFUN 
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INTEGER  ENDCON, ERROR 
LOGICAL  NOCONT, ERR, LEFT 
DIMENSION  TM  P  (  2  ) 

DIMENSION  CP (2, 200) 

COMMON /LINE/ BASE  (2  )  ,  DIR  (2  ),  DF  DS/ S  PEC  S  /  R  E  S  ,  DELTA  ,  TH  ETAN  ,  CON- 
COMMON /L  IMIT/YMIN , FMIN 

COMMON  /IOUNIT/  I  0 1 N , I OOUT , I S PT PF , LOUN I T , N D I R TU , NB SC AT 
DATA  FMIN/-29. / 

MAXDIM-200 
ERROR-O 
YM 1N-YMN 
C 

C  FINDING  THE  POINT  ON  THE  CONTOUR, OR  IF  A  CONTOUR  EVEN  EXISTS. 

C 

N  S  ER  C  H  =  0 

CALL  CL IMB(FUNCT,GFUN, CP, FCP.NSERCH, NOCONT) 

C  ***  APPROPRIATE  ACTION  IF  THE  CONTOUR  DOES  NOT  EXIST, 

IF (NOCONT  )GO  TO  99 
L-  1 

LEFT-.  TRUE. 

CALL  PTFIND (L E FT , C F UN , M AXD IM  ,  C P , L, ENDCON ) 

C 

C  **  THE  NEXT  TWO  DO  LOOPS  ARE  FOR  SWITCHING  THE  POSITIONS 
C  OF  THE  POINTS  IN  THE  ARRAY  AROUND  FOR  EASE  IN 
C  WORKING  WITH  LATER  ON. 

C 

MI D=L 12 
DO  17  I  *  1 ,  M  I D 
J  -L  -I  +  1 
DO  19  K ■ 1  ,  2 
TMP  (K)«CP(K,  I  ) 

C  P(K  ,  lKP(K.J) 

C  P  ( K  ,  J  )  -I  M  P  ( K  ) 

19  CONTINUE 
17  CONTINUE 

C  **  IS  THE  CONTOUR  CLOSED  OR  IS  THE  ARRAY  FILLED 
IF (ENDCON  )  9  7,  7  6,  2  7 
27  IF  (ENDCON.  EQ.  2)ERROR-2 

I F (ENDCON. EQ. 2 ) WRITE (I OOUT, 796 ) 

LEFT-. FALSE. 

CALL  PTFIND (LEFT ,G FUN, MAX DIM , CP, L, ENDCON ) 

MAXDIM-L 

IF (ENDCON )  9  7,  7  6,  95 
76  CONTINUE 
M  AXD  IM  -L  + 1 

CP(  1  ,  MAXDIM  )  -C  P(  1  ,  1  ) 

CF(2,MAXDIM)-CP(2,  1  ) 

GO  TO  96 
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95  IF(ENDCON.  EQ.  1  )C0  TO  98 

96  ERROR-2 
MAXD1M-L 

WRITE (IOOUT ,  796  ) 

796  F  ORMAT ( 5  3H  ***  DIRTRAN-I  ERROR  -  CONTOUR  TRACING  ROUTINE  ST'TK  / 

1  56H  CHECK  FOR  DISCONTINUITY  OF  FUNC7  V  ) 

GO  TO  999 

97  ERR0R--1 

WRITE (IOOUT, 797) 

797  FORMAT (5 1H  ***  DIRTRAN-I  ERROR  -  ARRAY  CPTS  NOT  LARCE  ENOUGH  ) 

GO  TO  999 

98  ERROR-O 
GO  TO  999 

99  ERROR® I 

WRITE (IOOUT,  7  9  9  ) 

7  99  F ORMAT ( 5 IH  ***  DIRTRAN-I  WARNING  -  CLOUD  CONTOUR  NOT  FOUND  / 

1  49H  MAY  HAVE  DISSIPATED  ) 

999  RETURN 
END 

***********************  SUBFILE  14  ***************************** 


FUNCTION  GF  UN  (S  ) 

IMPLICIT  INTEGERS  ( I -N  ) 

GF  UN  IS  THE  RESTRICTION  OF  THE  TWO  DIMENSIONAL  FUNCTION,  F,  TO 
A  LINE.  I.E.  FORM  C(S)-F(X,Y),  WHERE  (X,Y)»BASE45*DIF. . 

EXTERNAL  FL'NCT 

DIMENSION  P  (  2  ) 

COMMON /LINE/ BASE (2), DIR(2),DFDS/S PECS/RES, DELTA, THE TAN, CON 
CALL  VS UM (BASE, DIR, S,P) 

GFUN-F  UNCT  (P  ( 1  )  ,  P  (2  )  ) 

RETURN 

END 


**************************  SUBFILE  15  ************************ 

THIS  MODULE  IS  A  SUBROUTINE  THAT  FINDS  A  POINT  ON  A  CONTOUR 
BY  FINDING  THE  GRADIENT  VECTOR  AT  THAT  POINT  AND  MARCHING  ALONG 
IT  UNTIL  IT  FINDS  ITSELF  IN  A  REGION  GREATER  THAN  THE  CONTOUR  LEVEL. 
AT  WHICH  POINT  IT  MARCHES  HORIZONTALLY,  HALVING  THE  STEP  SIZE 
UNTIL  THE  CONTOUR  IS  REACHED  WITHIN  SPECIFIED  RESOLUTION. 

IN  ADDITION  IT  WILL  DETERMINE  IF  A  CONTOUR  EXISTS. 


C 

C  ARGUMENTS  PASSED. 
C 
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C  INPUT 

C  FUNCT-THE  FUNC T ION (X , Y  )  ALSO  GIVEN  IN  EXTERNAL. 

Pl-THE  STARTING  POINT. 


OUTPUT 

PI  -  THE  POINT  ON  THE  CONTOUR  OR  THE  POINT  AT  WHICH 
THE  FUNCTION  REACHES  A  MAXIMUM  BELOW  THE  CONTOUR 
LEVEL 

FP1  -  THE  VALUE  OF  THE  FUNCTION  AT  P 
NOCONT-THE  ERROR  FLAC. 

F-NO  PROBLEM 

T-NO  CONTOUR  FOUND. 

ERR-ERROR  FLAG  RETURNED  BY  'NTRSCT' 

F-NO  ERROR 

T-1TERATION  D1VERCED  OR  MAXIMUM  SEARCH  AREA  EXCEDF.D 
IN  ADDITION, IN  COMMON  ARE.., 

YMIN-THE  LOWER  LIMIT  ON  Y. 

DELTA-  THE  STEP  S I Z E , MO D I F I  ED  IN  THIS  SUBROUTINE. 
CON-THE  CONTOUR  LEVEL. 

RES-THE  RESOLUTION  LENGTH 

OTHER  VARIABLES  INCLUDE 

GRAD-THE  GRADIENT  VECTOR 

PO-THE  CURRENT  POINT  ON  THE  GRADIENT. 

Pl-THE  POINT  ON  THE  GRADIENT  BEING  TESTED 
TO  SEE  ABOUT  CONTOUR  EXISTENCE. 

FPO, FP1-THE  FUNCTION  VALUES  OF  PC  AND  PI. 

CALLED  SUBROUTINES 


C 

c 


GRAD2-FINDS  THE  GRADIENT  VECTOR  OF  A  FUNCTION  AT 
A  POINT  AND  THE  SLOPE  THERE. 

UNIT-CALCULATES  THE  NORM  AND  MAGNITUDE  OF  A  2  VECTOR. 
VSUM-VECTOR  SUM  OF  THE  FORM  C-A+SB  WHERE  S  IS  SCALAR 
MULTIPLIER  OF  B. 


SUBROUTINE  CLIMB(FUNCT,CFUN,PI,FPl,NSERCH,NCCONT) 

IMPLICIT  INTEGER*4  (1 -N  ) 

EXTERNAL  FUNCT, CFUN 
LOGICAL  NOCONT 

DIM EN SION  CRAD(2),P0(2),P1(2),P(2) 

COMMON /LINE/ BASE(2),DIR(2),DFDS/S PECS/RES, I  LTA  ,  THE TAN , CON 
C  OMMON /LIMIT/YMIN.FMIN 
NOCONT-. FALSE. 
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0NEM--1. 0 

IF  (NSERCH.  EQ.  0)G0  TO  7 
D  ELTA-S IGN (DELTA, FLOAT (NSERCH)  ) 

FP1«FUNCT(P  1  (1  )  ,  PI  (2  )  ) 

IF(FP1. LT. CON )CO  TO  25 
GO  TO  22 

3  CONTINUE 

PO (1 )-P  1  (1  ) 

P  0 ( 2 )-P  1  (2  ) 

FPO-F  P  1 

C  **  FINDING  THE  UNIT  GRADIENT  AND  THE  NEXT  POINT  ALONG  IT. 

4  CALL  GRAD2 (PO, FUNCT, RES , GRAD, DFDS  ) 

5  CALL  VS UM(P0, GRAD, DELTA, P  1  ) 

C  **  IS  THE  POINT  HEADING  BELOW  YMIN  ** 

IF  (P  1  (2  )  .  GE.  YMIN  )GO  TO  7 
P  1  (2  )  **YM  IN 

CALL  VSUM(P  1,  PO,  ON  EM, GRAD) 

CALL  UN  IT (GRAD , GRAD, DELTA ) 

1F(ABS(DELTA).LT.RES)G0  TO  25 

7  FP1«FUNCT(P  1  (1  )  ,  P  1  (2  )  ) 

C  **  HAS  THE  CONTOUR  BEEN  CROSSED  ** 

8  I  F (  F P  1 .  C E .  CON  )G 0  TO  22 
IF  (FP1.  GT.  FPO)GO  TO  3 
DELTA-DELTA/2. 

IF(ABS (DELTA ). LT. RES )GO  TO  25 
GO  TO  5 

25  NOCONT-. TRUE . 

GO  TO  99 
22  CONTINUE 

C  BEGIN  HORIZONTAL  SEARCH 
PO (2 ) -P  1  (2  ) 

31  P  0 ( I  ) -P  1  ( 1  ) 

FPO-F  P  1 

40  P 1 (1  )-PO(l  ) +D  ELTA 

FP  1-FUNCT  (P  1  (1  )  ,  PI  (2  )  ) 

1 F ( ABS (DELTA ) . LT. RES / 2 . )C0  TO  99 
IF  (FP1.  CE.  CON  )C0  TO  31 
DELTA-D  ELTA/2 . 

GO  TO  40 
99  CONTINUE 
RETUR  N 
END 

SUBFILE  16  ********************** 


THIS  SUBROUTINE  CALCULATES  THE  UNIT  GRADIENT  VECTOR 
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C  OF  A  FUNCTION  AT  THE  GIVEN  POINT  USING  THE  FORMULA 
C  PARTIAL  DF/DX  -  (F ( X 4fi , Y ) -F ( X , Y )  )  /R  WHERE  R  IS 
C  SMALL,  SIMILARLY  FOR  PARTIAL  DF/DY.  IT  THEN 
C  NORMALIZES  THE  RESULTANT  VECTOR  FOR  THE  UNIT  CRADIENT 
C  AND  FINDS  THE  SLOPE, WHICH  IS  THE  MAGNITUDE  OF 
C  OF  THE  REGULAR  CRADIENT. 

C 

ARGUMENTS  PASSED 

PT-THE  POINT  AT  WHICH  THE  UNIT  GRADIENT  IS  FOUND. 
FUNCT-THE  FU  NC  TI  0  N  (X  ,  Y  ) 

RES-R 

GRAD-THE  UNIT  GRADIENT 
SLOPE -SLOPE  AT  PT 

OTHER  VARIABLES 

COO, C 1 0, CO  1-THE  FUNCTION  AT  THE  POINTS  F(X,Y),F(X+R,Y), 
AND  F  (  X  ,  Y  +R  )  RESPECTIVELY. 

SUBROUTINES  CALLED 

UN  IT-NORMAL IZ ES  A  VECTOR  AND  FINDS  ITS  MAGNITUDE. 

SUBROUTINE  CR AD  2 ( P T , FU NC T , R E S , C R AD , S LO PE ) 

IMPLICIT  INTECER*4  (I-N) 

DIMENSION  PT  (2 ) , GR AD ( 2 ) 

COO-FUNCT (PT  ( 1  ) , PT (2  )  ) 

C  1  0-FUNC  T (PT  (  1  )+RES , PT (2 )  ) 

CO  1  -FUNCT  (PT  (1  )  ,  PT  (2  )+RES  ) 

GRAD  (1  )-  (C  1  0-C00)  /RES 
GRAD  (2  )  >=  (C01-C00)  /RES 
CALL  UNIT (GRAD, CRAD, SLOPE  ) 

RETURN 
END 


*******************************  SUBFILE  17  ******************** 

THIS  SUBROUTINE  FINDS  THE  INTERSECTION  OF  A  LINE  SPECIFIED  BY 
A  BASE  POINT,  BASE,  AND  A  DIRECTON  VECTOR,  DIR.  THE  ROUTINE 
DETER!'.  I N  E  S  A  REAL  NUMBER,  S,  SUCH  THAT  BASE+S*DIR  IS  THE 
POINT  OF  INTERSECTION  WITH  A  CONTOUR  THROUGH  N EW T ON-R A  PH SON 
ITERATION . 


C  INPUT  TERMS. 

C 

C  SLIM-THE  LIMITING  VALUE  OF  S,  TO  PREVENT  THE  INTERPOLATION 

C  FROM  BREAKINC  DOWN  IN  UNUSUAL  CASES  AND  COINC  INTO  AN 

C  INFINITE  LOOP 

C  GFUN -FUNCTION  OF  S  EQUIVALENT  TO  F U NC T I  ON (X , Y ) . 
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c 

OUTPUT  TERMS. 

C 

C  S-A  SCALAR  PARAMETER  SPECIFYING  POINTS  ON  THE  LINE. 

C  P-THE  POINT  OF  INTERSECTION. 

C  ERR-  AN  ERROR  FLAG. 

C  F-NO  PROBLEM 

C  T-NO  INTERSECTION  FOUND  WITHIN  THE  SET  LIMITATIONS. 

C 

C  TERMS  IN  COMMON. 

C  BASE-  A  2-VECTOR  CONTAINING  THE  COORDINATES  OF  THE  BASE 

C  POINT  OF  THE  LINE.  BASE  MUST  BE  INITIALIZED  BEFORE 

C  THE  CALL  TO  NTRSCT. 

C  DIR-  A  UNIT  DIRECTION  2-VECTOR  FOR  THE  LINE.  DIR  ALSO  MUST 

C  BE  INITIALIZED  PREVIOUS  TO  THE  CALL. 

C  CON-THE  CONTOUR  LEVEL  TO  BE  FOUND. 

C  DFDS-  DF/DS  WITH  F  «F  UNCT  10  N  (X  ,  Y  )  ,  USED  AS  AN  INITIALIZER 

C  FOR  THE  INTERPOLATION.  FIRST  ONE  MUST  BE  INITIALIZED 

C  BEFORE  THE  CALL. 

C  RES-THE  RESOLUTION 

l. 

C  CALLED  SUBROUTINES 
C 

C  VSUM-VECTOR  ADDITION  SUBROUTINE  OF  FORM  VS UM ( A , B  ,  S  ,  C )  *= 

C  C-A+S*B  WHERE  A,8,C  ARE  2-VECTORS  AND  S  SCALAR. 

C 

C 

SUBROUTINE  NTRSCT(S,P,CFUN,  SLIM, ERR) 

IMPLICIT  INTEGERS  (I-N) 

EXTERNAL  GFUN.FUNCT 
LOCICAL  ERR 
DIMENSION  r(  2) 

COMMON  /L INE/BASE(2),DIR(2),DFDS/S PECS/RES, DELTA, THETA N, CON 
COMMON /L IM IT/YMIN, FMIN 
SOLD-0. 0 
ERR=. FALS  E. 

C 

C  ***  ITER  IS  THE  NUMBER  OF  ITERATIONS.*** 

1TER-0 

F  SO  LD  *G  F  UN (SOLD) 

IF (FSOLP . LE. FMIN)GO  TO  999 
DS-  (CON-FSOLD) /DFDS 
S-SOLD+DS 

IF  (ABS(S).LT.SLIM)GO  TO  110 
DS-S1CN (SLIM/2.  , D  S ) 

S -SOLD+D  S 
1  1  0  I TER-I TER+1 


62 


o  r> 


FS-G  FUN (S  ) 

IF (FS. LE. FMI N )G  0  TO  999 
5  DFDS"(FS-FSOLD)/ (S  -S  OLD  ) 

D  S  *  (CON-FS  ) /DFDS 

SOLD«S 

FSOLD-F  S 

S-S+DS 

IF (ABS (DS ) . LT.RES )GO  TO  998 
IF (ABS (S ) .GT. SLIM)CO  TO  999 
I F ( ITER. GT. 9  )G  0  TO  9  99 
GO  TO  110 

9  98  CALL  VSL'M(BASE,DIR,S,P) 

CO  TO  1000 
999  ERR*. TRUE. 

1000  RETURN 
END 
C 

C  **************************  SUBFILE  18  ************************* 

C 

C  PTFIND  DETERMINES  A  STRING  OF  POINTS  ON  THE  CONTOUR  MARCHING 

C  EITHER  CLOCKWISE  (.NOT. LEFT)  OR  COUNTERCLOCKWISE  (LEFT)  UNTIL 

C  THE  CONTOUR  CLOSES,  A  BOUNDARY  IS  MET,  THE  ARRAY  CPTS  IS  FILLED, 

C  OK  AN  ERROR  OCCURS.  A  STARTING  POINT  MUST  BE  PROVIDED.  PTFIND 

C  BEGINS  BY  CALCULATING  THE  TANGENT  TO  THE  CONTOUR  WHICH  IS 
C  PERPENDICULAR  TO  THE  CR AD  I  ENT .  SUBSEQUENT  POINTS  ARE  FOUND  BY 

C  EXTENDING  THE  TANGENT  THROUCH  THE  TWO  MOST  RECENTLY  DETERMINED 

t  POINTS  AND  SEARCHING  IN  THE  P E RP E N D UC U LA R  DIRECTION.  IN  CASE 

C  OF  ERROR,  ONE  RESTART  IS  ALLOWED. 

C  THE  STEP  SIZE,  DELTA,  IS  HALVED  WHEN  THE  EXTRAPOLATED  POINT 

C  EXCEDES  BOUNDARY  OR  THE  SEARCH  IN  THE  PERPENDICULAR  DIRECTION 

C  DID  NOT  YIELD  A  POINT  WITHIN  LIMITS  DETERMINED  BY  CURVATURE. 

C  DELTA  IS  LENGTHENED  INVERSELY  AS  THE  LOCAL  CURVATURE  OF  THE  CURVE. 

C 

C  INPUT  TERMS 

C  GFUN-THE  ONE  VARIABLE  EQUIVALENT  TO  THE  FUNCTION 

C  RESTRICTED  TO  A  LINE  THROUGH  POINT,  BASE,  IN  D I R EOT  ION , D I R 

C  MAXDIM-THE  MAXIMUM  NUMBER  OF  POINTS  IN  THE  ARRAY. 

C  CP-THE  ARRAY  OF  CONTOUR  POINTS  AS  IT  INITIALLY  IS. 

C  L-THE  INDEX  OF  HOW  MANY  POINTS  HAD  BEEN  FOUND  BEFORE. 

C 

C 

C  OUTPUT  TERMS 

C  CP-CURRENT  ARRAY  OF  CONTOUR  POINTS 

C  L -CURRENT  INDEX  OF  POINTS.  ON  RETURN,  L  IS 

C  THE  NUMBER  OF  POINTS  FOUND. 

ENDCON-A  NUMERICAL  FLAC  THAT  TELLS 

WHY  THE  SUBROUTINE  IS  ENDING  THE  SEARCH. 

C  ENDCON  EQUALS. 
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-1«>THE  ARRAY  CP  IS  FILLED. 

0 ■ >T  H  E  CONTOUR  IS  CLOSED. 

I->THE  CONTOUR  HAS  MET  A  BOUNDARY. 

2  *>R  EQU  I R  ED  STEP  SIZE  BECAME  LESS  THAN  RESOLUTION- 

OTHER  PERTINENT  VARIABLES  IN  COMMON... 

CON-THE  CONTOUR  LEVEL 

BASE-THE  POINT  ALONG  THE  TANGENT. 

DFDS-T  HE  SLOPE  OF  THE  CURVE.  MODIFIED  IN  'NTRSCT' 

RES-THE  RESOLUTION  LENGTH 

DELTA-THE  STEP  SIZE  ALONG  THE  TANGENT, MODI F 1  ED  HERE. 
THETAN-THE  TANGENT  OF  THETA,  THE  ANGLESPREAD  OF  PERMISSIBLE 
SEARCH  AREA.  THE  SMALLER  THETAN  IS,  THE  SMOOTHER  THE 
POLYGON  OF  CALCULATED  POINTS  WILL  BE. 

YM IN-LOW  ER  BOUND  FOR  SECOND  COORDINATE 

OTHER  VARIABLES 

PI-CURRENT  POINT  ON  CONTOUR 

PNEXT-NEXT  POINT  ON  CONTOUR, NOT  YET  APPROVED  BY  'ENDTST'. 
FULTAN-THE  APPROXIMATION  OF  THE  TANGENT  AT  PNEXT  BY 
THE  VECTOR  FROM  PI  TO  PNEXT. 

TAN-THE  UNIT  VECTOR  OF  FULTAN 

DIST-THE  MAGNITUDE  OF  FULTAN, USED  IN  'ENDTST'. 

S-THE  DISTANCE  ALONG  THE  GRADIENT  FROM  BASE  TO 
A  POINT  ON  THE  CONTOUR. 

SMAX-THE  LARGEST  ALLOWABLE  VALUE  OF  S. 


SUBROUTINES  CALLED. 

VS UM ( A , B , S , C ) -VECTOR  ADDITION  OF  THE  FORM  C«A+S*B  WHERE 
A , B , C  ARE  2-VECTORS  AND  S  IS  SCALAR. 

PERP(A,B)-A  ,B  ARE  2-VECTORS  AND  B  IS  A  VECTOR  ROTATED 
90  DEGREES  COUNTERCLOCKWISE. 

NTRSCT  (S  , P, CFUN , SLIM , I C , ERR  )  -F  INDS  THE  INTERSECTION  OF  A 
LINE  DRAWN  FROM  BASE  (IN  COMMON)  ALONG  DIR, THE  DIRECTION 
VECTOR  (IN  COMMON),  WITH  THE  CONTOUR  LEVEL  (IN  COMMON). 

S  IS  THE  DISTANCE  FROM  BASE  TO  THE  CONTOUR, P  IS  THE 
C  POINT  OF  INTERSECTION,  CFUN  IS  CFUN(S)  AND  IS  A  ONE 

C  VARIABLE  EQUIVALENT  OF  THE  FUNCTION  WHOSE  CONTOUR 

C  IS  BEING  SOUGHT,  SLIM  IS  THE  MAXIMUM  ALLOWABLE 

C  VALUE  OF  S,  AND  ERR  IS  THE  ERROR  FLAG (L OC I C AL  )  . 

C  IC  IS  A  LOGICAL  VARIABLE  WHERE  T  MEANS  FOR 

C  NTRSCT  TO  STOP  AFTER  5  ITERATIONS. 

C  ENDTST-F (P 0, P, TAN , D1 ST , SMAX, I , M, ENDCON ) -F I NDS  IF 

C  THE  CONTOUR  IS  CLOSED,  HAS  RUN  INTO  THE  Y  BOUNDARY 

C  OR  IF  THE  ARRAY  OF  CONTOUR  POINTS  IS  FILLED. 
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C  UNIT (A , B , S ) -A , B  ARE  2-VECTORS  AND  B  IS  THE  NNORM  OF  A 

C  AND  S  IS  THE  MAGNITUDE  OF  A. 

C 

C 

C 

SUBROUTINE  PTFIND (L E FT , G F UN , M AXD IM , CP , L , ENDC ON ) 

IMPLICIT  I  NTECER*4  (I-N) 

EXTERNAL  GFUN, FUNCT 
INTEGER  ENDCON 
LOGICAL  ERR, ONCE, LEFT 
DIMENSION  CHORD  (2  ),  TAN  (2  ) 

DIMENSION  C  P  (  2  ,  MAXD  IM  ) 

COMMON/L INE/BASE (2), DIR(2),DFDS/S  PECS/RES, DELTA, THE TAN, CON 
COMMON /L IM IT/YMIN , FMIN 
ERR- . FALSE. 

RES2-RES/2. 

C 

C  INITIALIZE  SEARCH  FOR  CONTOUR  POINTS  MAKINC  USE  OF  THE  FACT 

C  THAT  THE  TANGENT  TO  THE  CONTOUR  IS  PERPENDICULAR  TO  THE 

C  GRADIENT  OF  THE  FUNCTION. 

C 

9  ONCE-. TRUE. 

CALL  GRAD2 (C  P( 1 , L) , FUNCT, 5. *RES , DIR , DFDS ) 

CALL  PERP (D IR , TAN ) 

I F ( . NOT. LEFT )GO  TO  31 
TAN  ( 1  )  --TAN  ( 1  ) 

TAN  (2  )  --TAN  (2  ) 

DFDS  «-DFns 
31  DELTA-10. 

TANTH- 1 . 

C  ***  STEPPING  DELTA  ALONG  THE  TANGENT  *** 

10  CALL  VSUM (CP(1 , L) , TAN, DELTA, BASE ) 

C  ***  IS  BASE  BELOW  THE  MINIMUM  PERMISSIBLE  Y*** 

1 F( BASE (2 ) . CT. YM IN)GO  TO  29 
DELTA-DELTA/2. 

IF (ABS (DELTA ). LT. RES2 )C0  TO  91 
GO  TO  10 

DETERMINE  THE  RANGE  OF  SEARCH  PERPENDICULAR  TO  THE  EXTENDED 
TANGENT,  SMAX,  AND  CALL  NTRSCT  TO  LOCATE  THE  CONTOUR  POINT 
C 

29  SMAX -AM  IN  1 (ABS (D ELTA*T ANTH+R  ES ) , ABS(  (BASE (2 ) -YMIN)  /D I R (2 )  ) ) 
DFDSO-DFDS 

30  CALL  NTRSCT (S , CP(1 , L+l ), CFUN, SMAX, ERR ) 

C  ***  HAS  THE  CONTOUR  BEEN  FOUND  WITHIN  LIMITING  CONDITIONS** 

I F ( . NOT. ER  R ) C  0  TO  3  2 

DFDS-DFDSO 

DELTA-DELTA/2. 
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C  CHECK  TO  SEE  IF  STEP  SIZE,  DELTA,  IS  BELOW  MINIMUM,  RES2. 

C  IF  SO,  ALLOW  ONE  RESTART  THEN  ABORT 

C 

IF(ABS (DELTA) .CE. RES 2)CO  TO  10 
IF (ONCE  )G0  TO  92 
GOTO  9 

C  DETERMINE  TANGENT  FOR  NEWLY  FOUND  POINT 

32  CALL  VSUM  (  C  P  ( 1  ,  L  +  l  )  ,  CP(1  ,  L)  ,  -  1.  ,TAN  ) 

CALL  UNIT (TAN , TAN, DIST) 

CALL  PERP (TAN , DIR ) 

C  **  IS  THE  CONTOUR  CLOSED  ** 

C  WHAT  THIS  NEXT  CHECK  DOES  IS  TO  LOOK  WITHIN  A  BOX 

C  DELTA  BY  2*SMAX,  WHERE  P  AND  THE  PREVIOUS  POINT  ARE  AT  EACH 

C  END  OF  THE  LINE  SEGMENT  BETWEEN  THE  MIDPOINTS  OF  THE  TWO 
C  OPPOSITE  SIDES  2*SMAX  LONG,  TO  SEE  IF  THE  STARTING  POINT 
C  IS  WITHIN. 

CALL  VSUM (CP, CP(1 , L) , -1.  .CHORD  ) 

VAR1 -DOT  PR  D(DIR, CHORD) 

V  AR2 -DOT  PR  D(TAN,  CHORD) 

IF ( ABS (VAR  1  ) . CE. RES*2 .  )G0  TO  35 
I  F  (  (VAR2.  LT.  DI  ST  )  .  AND.  (  VAR2.  GT.  0.  0  )  )C  0  TO  90 
C  BEGIN  SEARCH  FOR  NEXT  POINT 

3  5  L-L  +  l 

IF(CP(2, L) . LT. RES 2+YM IN )CO  TO  91 
IF  (L  .  GE.MAXD  IH  )CG  TO  69 

DELTA-AM1N  1(1  0.  .DELTA*  (SMAX  +  2.*RES)  /  (  (ABS(S  )+RES)  *2.  )  .DELTA *2. 
ONCE-. FALSE. 

TANTH-THETAN 
GOTO  10 

89  ENDCON--1 
GO  TO  99 

90  ENDCON-O 
L-L  +  l 
GO  TO  99 

91  ENDCON-i 
GO  TO  99 

9  2  E  N  DC  ON  -  2 
9  9  RETURN 
END 

*******************„****  SUBFILE  19  *************************** 

C 

SUBROUTINE  VS UM ( A , B , S , C ) 

IMPLICIT  INTECER*4  (I-N) 

D IMEN  SION  A(2),B(2),C(2) 

C  ***  C-A+S*B  WHERE  A,B,C  ARE  VECTORS  AND  S  IS  SCALAR 
DO  14  J- 1 , 2 
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14  C  ( J  )-A ( J )+S  *B ( J  ) 
R  ETURN 
END 


SUBFILE  20 


SUBROUTINE  PERP(A,B) 

IMPLICIT  INTECER*4  (I-N) 

DIMENSION  A ( 2  )  , B ( 2  ) 

C  ***  B  IS  ROTATED  90  DECREES  COUNTERCLOCKWISE  FROM  A 
B  ( 1  )— A  (2  ) 

B  (  2  )  -A  (  1  ) 

R  ETUR  N 
END 
C 

C  *************************  SUBFILE  21  *************************** 

C 

SUBROUTINE  U N 1 T  ( A , B , XN ORM ) 

IMPLICIT  INTECER*4  (I-N) 

DIMENSION  A  (2  )  ,B(2  ) 

C  ***  B  IS  THE  NORM  OF  A,  AND  XN  ORM  IS  THE  MAGNITUDE 
XNORM*SQRT(A(l  )**2+A(2)**2) 

B ( 1  ) -A  ( 1  )  /XN  ORM 
B  (  2  )  -A  (  2  )  /XN  ORM 
RETURN 
END 
C 

C  **************************  SUBFILE  22  ************************ 

C 

FUNCTION  DOTPRD(A.B) 

IMPLICIT  INTEGER *4  (I-N) 

DIMENSION  A ( 2  )  ,  B  (  2  ) 

C 

C  DOT  PRD  IS  THE  SCALER  PRODUCT  OF  A  AND  B 

C 

DOT  PR  D*=A  (  1  )  *B  (  1  )  +A  (  2  )  *B  (  2  ) 

R  ETURN 
END 
C 

C  *************************  SUBFILE  23  ************************ 

C 

C  THIS  SUBROUTINE  TAKES  AN  ARRAY  OF  POINTS  OF  A  CONTOUR 

C  AND  FINDS  THE  HIGHEST  POINT,  THE  TWO  POINTS  ALONG  A 

C  LINE  OF  SICHT  (THESE  TWO  MAY  NOT  BE  IN  THE  ARRAY  AND  WILL 

C  BE  FOUND  BY  INTERPOLATION  FROM  RATIOS)  AND  THEIR  CENTROID, 

C  THE  CENTROID  AND  THE  TWO  POINTS  WITH  THE  SAME  Y  VALUE,  THE  LENGT 
C  OF  THE  LINE  OF  SIGHT,  AND  THE  SHEAR  DISTANCE.  NOTE-  THE  POINT  ON 
C  THE  CENTROID  LINE  WITH  THE  GREATER  X  VALUE  IS  THE  LEADINC 
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C 
C 

C  INPUT  TERMS 
C 

C  SPCHT-A  SPECIFIED  HEIGHT  AT  WHICH  THE  WIDTH  IS  DESIRED 

C  CP-THE  ARRAY  OF  CONTOUR  POINTS 

C  L-THE  NUMBER  OF  POINTS  IN  CP. 

C 

C 

C  OUTPUT  TERMS 
C 

C  Z II  I  -  Z  COORDINATE  OF  HIGHEST  POINT  ON  CONTOUR. 

C  XHI-  X  COORDINATE  OF  HIGHEST  POINT  ON  CONTOUR. 

C  XCL-  X  COORDINATE  OF  LEFTMOST  POINT  ON  CONTOUR. 

C  XCR-  X  COORDINATE  OF  RIGHTMOST  POINT  ON  CONTOUR. 

C  CNTRD-T HE  CENTROID 

C  CENWTH-  WIDTH  AT  THE  CENTROID  HEIGHT 

C  XS  PC  L-  X  COORDINATE  OF  LEFTMOST  POINT  OF  CONTOUR  AT  THE 

C.  SPECIFIED  HEIGHT,  SPCHT 

c  XS  PC  R-  X  COORDINATE  OF  RIGHTMOST  POINT  OF  CONTOUR  AT  THE 

L  SPECIFIED  HEIGHT 

C  SPCWTH-  THE  WIDTH  OF  THE  CONTOUR  AT  THE  SPECIFIED  HEIGHT 

C 

C  SUBROUTINES  CALLED 
t 

C  I NTRP-PERFORMS  A  SIMPLE  INTERPOLATION  BY  RATIO  TO  FIND 

C  A  POINT  WITH  A  GIVEN  Y  COORDINATE,  GIVEN  A  POINT 

C  ON  EITHER  SIDE.  THE  CALL  IS 

C  CALL  INTRP  (P  1,  P2,  Y,  P)  WHERE  P1.P2  ARE  THE  GIVEN 

C  POINTS,  Y  THE  GIVEN  Y  COORD.,  AND  P  THE  POINT  RETURNED. 

C 

C  IN  INTP.P...  VSUM  (  A  ,  B  ,  S  ,  C)  WHERE  A,B,S  ARE  GIVEN  AND 

C  C*A+S*B.  A , B , C  ARE  2-VECTORS  AND  S,  A  SCALAR. 

C 

C 

SUE  ROUT  1  NE  PT  S  PE  C  (L  ,  C  P ,  XH  I  ,  ZH  I  ,  XC  L  ,  CNTRD,  XC  R  ,  C  ENW  T  H  ) 

IMPLICIT  INTEGER*^  (I-N) 

DIMENSION  CNTRD (2 ) , CP (2 , L ) 

LOGICAL  ERR.BTWN 

BTWN  (A,E,C ) “  (A . LE.  B. AND.  B.LT.  C)  .OR.  (A.CE.B.AND.  B.CT.  C) 

ERR*. FALSE. 

N  H  I  =  1 
N  LO*  1 
N  L  *  1 
N  R  =  1 

DO  10  J-  1  ,  L 

IF(CP(2,J).CT.CP(2,NHI))NHI-J 
I  F  (CP(2  ,  J  )  .  LT.  CP  (2  ,  NLO  )  )N  LO-J 
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1F(CP(1,J).LT.CP(1,NL))NL-J 
IF(CP(1,J).CT.CP(1,NR))NR-J 
10  CONTINUE 

XHI-CP(1  ,  N II 1  ) 

ZH1-CP(2,NH1  ) 

CNTRD(2)«A!1AX1(CP(2,NL),CP(2,NR)) 

NCR  *N  L 
NCL«NR 
DO  2  0  J  «  1  ,  L 

IF  (  .  NOT.  BTWN  (CP(2  ,  J  )  ,  CNTRD  (2  )  ,  CP(2  ,  J  +  l  )  )  )G  0  TO  20 
IF(CP(1,J).LT.CP(1,NCL))NCL«J 
I F(C P ( 1  , J ) . CT. CP( 1 , NCR)  )NCR-J 
2  0  C  0  N  T 1 N U E 

CALL  1NTRP (CP(1 , NCL) ,  C  P  ( 1 , NCL+1 ),CNTRD(2 ) ,XCL) 

CALL  INTRP(CP( 1 , NCR) ,CP(1 , NCR+1 ) ,CNTRD(2 ) ,XCR) 

CNTRD  (  1  )  ■=  (XC  L  +  XC  R  )  /2  . 

C  ENWTH  »XC R-XC  L 

C  I  F  (  .  NOT  .  BTWN  (2LO  ,  SPCHT,  ZH  I  )  )CO  TO  99 

c  ncr  =:;  l 

C  NCl.-NK 

C  DO  3  0  J  *=  1  ,  L 

i  IF (. NOT. BTWN (CP(2 , J ), SPCHT, CP(2, J  +  l ))  )GO  TO  30 

C  IF(CP(1,J).LT.CP(1,NCL))NCL-J 

C  1 F ( C  P( 1 , J ) . CT. CP( 1 , NCR)  )N  CR«J 

t  J  0  C  C  N  T  I  N  L  E 

t  CALL  1 NTRP (CP( 1 , NCL) ,CP( 1 , NCL+ 1 ), SPCHT, XSPCL) 

C  CALL  1NTRP (CP( 1 , NCR ) ,CP( 1 , NCR+ 1) , SPCHT, XS PCR) 

C  S TCWTH»XS PCR-XS PC L 

C  C  C  T  0  9  9  9 

C  9  l<  ERR=.  TKl'l  . 

(  XS  IV  L-D.  0 

t  X  S  PC  R  =  0  .  0 

i  S  PC  W1  II  =  0  .  0 

9  9  R  r.  1  L  R  X 
E  N  I' 

L 

(  ***********************  Sl'PFILE  2  U  ************************** 

SOB  hoc:  1 N I  1 N  T  R  P  (  P  1 ,  P  2 , V , X  ) 

IMPLICIT  I  NT  EGER *U  ( I  -N  ) 

(  *  ****************************************************** 

c 

t  GIVEN  POINTS  PI  AND  P2  AND  SECOND  COORDINATE,  Y, 

c:  INI  R  P  DETERMINES  THE  FIRST  COORDINATE,  X,  OF  A  POINT,  ( X ,  V  I  , 

l  WHICH  IS  ON  A  LINE  DRAWN  THROUGH  PI  AND  P2 

( 

t  ************************************************************ 

DIMENSION  P 1  ( 2  )  , P  2 ( 2  )  , P ( 2  )  , D  1  F ( 2  ) 
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1  F  (  ABS  (P  2  (2  ) -P  1  (2  )  ).  LT.  1 .  E-30)C0  TO  89 
RATIO  -  (Y-Pl  (2  )  )  /  (P  2  ( 2  )-P  1  (2  )  ) 

0  N  EM - - 1 . 

CALL  VSUM(P2,  Pl.ONEM.DIF) 

CALL  VS UM (P 1, DIF, RATIO, P) 

X-P  (  1  ) 

CO  TO  99 

8  9  X-P  1(1) 

9  9  R  ETUR  N 

END 

CALCULATION  OF  O-ORDER  AND  INTERPOLATION  OF  HIGHER  ORDER  MOMENTS 

*********************  SUBFILE  25  ******************************* 


SUBROUTINE  MOM  ENT (V CRA V, Z IN , H , T I N , Q , XB AR , S  I  CW 2 ,  S 1 C P 2 ) 

IMPLICIT  I NTEGER*4  (I-N) 

R  EAL  M , N , NM 

DIMENSION  AL (9 ) , Z( 9 ) ,T (9 ) , XB ( 8 1 , 4, 4 ) , SU( 8 1 , 4 , 4 ) , SP ( 8  1  , 4, 4 ) , NM( 9  ) 
DIMENSION  VA  L ( 1 6 ) , X  VAL ( 8 ) , W ( 8 ) , X  I ( 4 ) , 1 B ( 4 ) , NTC ( 4  )  , I  I  (  4  )  , X  (  9 ,  4  ) 
LOGICAL  FIRST 

common  /wndprm/dxzo, dyxo, dzo, uo, m, n, zinv 


COMMON  /I  OUNIT/  I  0  I  N  ,  I  OOL'T,  I  S  PT  PF  ,  LOL’NI  T,  N  D  1  R  TU  ,  N  B  SC  AT 
EQU I VA LENC E  (Z(l  )  , X  ( 1 ,  1)  ),(T  (1  )  , X  (  1 , 2  )  )  ,  ( A  L  (  1  )  , X  (  1 ,  3  )  ) 
EQUIVALENCE  (NM ( 1  )  , X  (  1 , 4  )  ) 


DATA  FI  R ST/. TR U E. /, IB/4, 4,1,  3/,  ITC/11/, HREF/1./ 

**************************************************************** 


PUR  POSE 


TO  CONVERT  PARAMETERS  TO  N 0 N D I M E N S I  0 N A L  FORM  AND  THEN  COMPUTE 
THE  ZERO  ORDER  MOMENT  AND  INTERPOLATE  FROM  TABULATED  VALUES  (  F 
THE  HIGHER  ORDER  MOMENTS 

INPUT 

VCRAV  THE  GRAVITATIONAL  SETTLING  VELOCITIES  OF  THE  PARTICLE 
IN  METERS  /  SEC 

ZIN  THE  HEIGHT  (METERS)  AT  WHICH  THE  MOMENTS  ARE  DES1RFD 

H  THE  HEIGHT  OF  RELEASE  OF  THE  PARTICLES  IN  METERS 

TIN  THE  TIME  IN  SECONDS  AFTER  RELEASE 

OUTPUT 

Q  THE  VERTICAL  CONCENTRATION  OF  PARTICLES  IN  PARTS/METER 

AT  HEIGHT  Z 

THE  DISPLACEMENT  (METERS)  IN  THE  X  (IE  WIND>  DIRECTION 
OF  THE  CENTER  OF  MASS  OF  PARTICLES  AT  HEIGHT  Z 


XB  AR 


n  o 


S  I  GW  2 


C 

c: 

c 

c 


S  IGP  2 


THE  SQUARE  OF  THE  STANDARD  DEVIATION  OF  THE  WINDWARD 
DISPLACEMENT  OF  THE  PARTICLES  AT  HEIGHT  Z  IN  MFTF.FS**2 
THE  SQUARE  OF  THE  STANDARD  DEVIATION  OF  THE  CROSS-WIND 
DISPLACEMENT  OF  THE  PARTICLES  AT  HEIGHT  Z  IN  METERS** 2 


I 


C 

c 

c 

c 

c 

c 

L 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c: 

c 

c 

c 

c: 

c 

c: 

c 

c 

c 

c 

c 


SUBROUTINES  CALLED 

DTERPS  PUTS  THE  NEEDED  VALUES  OF  THE  TABULATED  MOMENTS 
INTO  A  ONE  DIMENSIONAL  ARRAY 

DTERri  A  FUNCTION  WHICH  RETURNS  THE  INTERPOLATED  VALUE 
FOR  GIVEN  ARGUMENTS  AND  ARRAYS 
CREEN  CALCULATES  THE  GREENS  FUNCTION  WHICH  IS  THE 
O-ORDER  MOMENT 


CALLED  BY  CW1ND 


******************: 


************************* 


I F ( . NOT . FIRST )GO  TO  5 


READ  IN  THE  TABLE  OF  MOMENTS  ON  THE  FIRST  CALL  OF  MOMENT 


Z  LOG  OF  NON-DIMENSIONAL  HEICHTS  AT  WHICH  MOMENTS  ARE  TABUIATE 

T  LOC  OF  NON-DIMENSIONAL  TIMES  AT  WHICH  MOMENTS  ARE  TAEl'LATF 

AL  NON-DIMENSIONAL  SETTLING  VELOCITIES  AT  WHICH  MOMENTS  AFF 

TABULATED 

NM  DIFFUS IVITY  POWER  LAW  EXPONENTS  AT  WHICH  MOMENTS  ARE 

TABULATED 

XU  TABULATED  VALUES  OF  LOGS  OF  FIRST  ORDER  MOMENTS  (RELATED 

TO  MEAN  HORIZONTAL  DISPLACEMENT) 

SW  TABULATED  VALUES  OF  LOGS  OF  WIND  SHEAR  COMPONENT  OF  SECOND 

ORDER  MOMENT  (CONTRIBUTES  TO  VARIANCE  IN  WIND  DIRECTION  • 

SP  TABULATED  VALUES  OF  LOCS  OF  SECOND  ORDER  MOMENT  COMMON  TO 

WIND  AND  CROSS-WIND  VARIANCES 


READ (N DIRT U, 1 ) 

1  FORMAT (4  I  3) 

NTC  (  I  )  -N  Z  -  I 
NTC (2  ) -N T-  1 
N  T  C  (  3  )  -N  A  -  1 
N  TC  (  4  )  -N  N  -  1 

R  EAD (N  DIRTU,  2  ) 

2  FORMAT (6E  1  3.  5  ) 
READ (NDIRTU, 2) 
R  EAD  (N  D  I  R  TU  ,  2  ) 
R  EA  D (N  D I R  TU , 2  ) 
NZT-N  Z*NT 


NZ ,NT,NA,NN 


(Z (  I  ) , 1- 1  ,  NZ) 

(T (I  ) , 1-1,  NT) 
(AL(I  )  ,  I-  1, NA) 
(N  M  (  I  ),I-1,NN) 
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DO  3  L-l, NN 

READ  (N  DIRT  U,  2  )  (  (XB  ( I  J,  K ,  L  )  ,  I  J  - 1 ,  NZT  )  ,  K  -  1  ,  NA  ) 

R  EA  D (N  D I RTU ,2)  ((SW(IJ,K,L),IJ-1,NZT),K«1,NA) 

READ  (NDIRTU,  2  )  (<SP(IJ,K,L),IJ-1,NZT),K-1,NA) 

3  CONTINUE 

FIRST-. FALSE. 

REWIND  NDIRTU 
5  CONTINUE 
C 

C  CONVERT  INPUT  PARAMETERS  TO  NO NDIM EN S 10 NAL  FORM 

C 

SCLU-DZO*H  ** (N -1 .  ) 

XI (1 )-Z IN/H 
XI  (2  )  -S  C  LU  *T  I  N  /  H 
XI (3 ) -VGRAV/SCLU 
XI  (4 ) -N 

CALL  GREEN(XI (1 ) , HREF, XI (2 ),XI (3 ),Q, IER) 

Q-Q/H 

IF (Q  . LE.  1. E-20)  GO  TO  999 

C 

C  TAKE  LOGS  FOR  LOGARITHMIC  INTERPOLATION 

C 

XI  (I  ) -A  LOG  (XI  (1  )  ) 

XI (2 ) -ALOG (X  I  (2  )  ) 

C 

C  DETERMINE  INDICES  OF  LOWEST  CORNER  POINT  OF  THE  CUBE  TO 

C  BE  USED  IN  INTERPOLATION  MAKINC  SURE  THAT  ENOUGH  CORNER  POINT 

C  OF  THE  CUBE  HAVE  TABULATED  VALUES 

t 

DO  100  1-1,4 
I  KI  ) -I  B  (  I  ) 

IOC'  CONTINUE 

DO  101  111-1,4 
1-5-1  I  1 

C  I  A  - 1  1(1) 

IF(XI(I)  .CE.  X  ( I A ,  I  )  .AND.  X(IA,I)  .LE.  X(IA+1,I))  C  ('  TO  101 

I  F  ( X  I  (1  )  .LT.  X  ( I  A  ,  I  )  .AND.  IA  .  EQ.  1)  CO  TO  101 

I F ( X  I  ( I  )  .CT.  X  ( 1  A  ,  I  )  .AND.  1A  .EQ.  NTC(I))  CO  TO  101 

I  SAV-I  1  (1  ) 

1  I  ( I  ) - 1  A  +  1  FI X(S ICN  (  1 .  ,XI  (I  ) -X (I  A , I )  ) ) 

IT-0 

DO  102  JI -1, 2 
JI X-J  1  +  11(1)  -  1 

DO  102  I J- 1 , 2 

I  JX-J  IX  +  (IJ  +  11(2)  -  2  )  *N  Z 
DO  102  K-  1  , 2 
KX-K -1  +  11(3) 

DO  102  L -  1  ,  2 
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LX-L-1  +  11(4) 

IF (XB (1 JX, KX,  LX)  .GT.  -100.)  IT-IT+1 

102  CONTINUE 

IF  ( I T  .GT.  ITC)  GO  TO  6 
I  I  (I  )-l SAV 
101  CONTINUE 
C 

C  PERFORM  THE  INTERPOLATION  WITH  DETERMINED  CUBE  OF  POINTS 

C 

DO  103  1-1,4 
12-1*2 
I  1-1  2-1 
I  A  - 1  1(1) 

XVAL  (I  1  )  -X  (1  A  ,  I  ) 

XVAL(I  2  )  -X  ( I  A+  1  ,  I  ) 

103  CONTINUE 

CALL  DTERPS  (I  I ,XB .VAL.NZ) 

XBAR-DTERPI (4, XI ,XVAL,VAL,-100.  ,W) 

CALL  DTERPS (I  I , SW, VAL, NZ) 

SIGW2-DTERPI (-4, XI , XVAL, VAL, -100. ,W) 

CALL  DTERPS  (I  I , SP, VAL, NZ) 

S  IGP2-DTERPI  (-4, XI , XVAL, VAL, -100.  , W) 

C 

C  CONVERT  THE  LOG  OF  THE  NONDIM ENS IONAL  VALUES  INTERPOLATED 

C  TO  THE  USUAL  DIMENSIONAL  FORM 

C 

SCL-U  0*H  ** (K+l.  )  /SCLU 
XB AR-SCL*E  XP (XB AR ) 

S  1GW2-SCL*SCL*EXP(S I  GW  2 ) 

S  IGP2-2.  *DXZ0*H*H  *EXP  (S  IGP2  ) 

S  I  GW  2  -S  IGW2+S  ICP2 
S  IGP2-D YX0*S 1 G  P  2 
9  99  RETURN 
END 
C 

q  ***********************  SUBFILE  26  ************************* 


C 

SUBROUTINE  DT E R P S ( I  I , X , VA L , NZ ) 

IMPLICIT  INTEGER*^  (I-N) 

DIMENSION  X ( 8 1 , 4 , 4 ) , VA  L ( 1  )  , I  I  ( 1  ) 

( 2  ************************************************************** 

C 

C 

C  PURPOSE 

C 

C  TO  SET  UP  A  ONE  DIMENSIONAL  ARRAY  OF  THE  VALUES  CORRESPONDING 

C  TO  THE  CORNERS  OF  THE  CUBE  WITHIN  A  TABULATED  ARRAY  WITH 

C  LOWEST  CORNER  INDICES  GIVEN 
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INPUT 


C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


II  SINGLY  DIMENSIONED  ARRAY  CONTAINING  THE  INDICES  OF  THE 

LOWEST  CORNER  OF  THE  CUBE 

X  A  TRIPLY  DIMENSIONED  ARRAY  CONTAINING  THE  TABULATED 

V“2  0 1ALUES  TO  BE  SET  UP.  THE  FIRST  INDEX  IS  THE  COLLAPSE! 
INDEX  FOR  THE  FIRST  TWO  INDICES  OF  A  FOUR  DIMENSIONAL 
ARRAY 

NZ  THE  RANGE  OF  THE  FIRST  INDEX  OF  THE  FOUR  DIMENSIONAL 

ARRAY 


OUTPUT 


C 

C 

t 

C 

C 

C 


VAL  SINGLY  DIMENSIONED  ARRAY  CONTAINING  THE  VALUES  OF  X 
FOR  THE  16  CORNER  POINTS  OF  THE  CUBE 


CALLED  BY  MOMENT 


A  A  A  A  A  i 


rAAAAAAAAAAAA 


11  =  0 

DO  10  4  L  =  1 ,  2 
L  X=L  +  11(4)  -  1 
DO  103  K  =  1  ,  2 
K»K  +  11(3)  -  1 


101 
102 
10  3 
1  04 


DO  102  Jl  =  1, 2 

J1X=  (J  1  +  11(2)  -  2  )  *NZ 

DO  101  1 J-l,  2 

I  JX=J  IX  +  IJ  +  11(1)  -1 

H-M  +  l 

VAL (M) -X (I JX, KX, LX) 

CONTINUE 

CONTINUE 

CONTINUE 

CONTINUE 

RETURN 

END 


C. 

C 

c 


C 

t; 

c 

c 

c 


**********************  SUBFILE  27 


:******* 


FUNCTION  DTERP1(NDIM,XI,XVAL,VAL,VMIN,W0RK) 
IMPLICIT  INTEGER*4  (I-N) 


************ 


*************** 


PUR  POSE 
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o  o 


C  PERFORMS  AN  N-D IMENSIONAL  LINEAR  INTERPOLATION 

C 

C 

C  INPUT 

C 

C  ND1M  -  THE  NUMBER  OF  DIMENSIONS.  (-  DON’T  RECALCUALTE  WEIGHTS) 

C  XI  -  THE  POINT  IN  THE  HYPERSPACE  AT  WHICH  THE  INTERPOLATED 

C  VALUE  IS  DESIRED.  XI  MUST  BE  A  VECTOR  OF  ATLEAST  N  D I  It 

C  IN  LENGTH. 

C  XVAL  -  THE  COORDINATE  VALUES  AT  THE  CORNERS  OF  THE  HYPERCUBE. 

C  THE  VECTOR  MUST  BE  SET  UP  LIKE  A  TWO-DIMENSIONAL  ARRAY 

C  (2  X  NDIM),  WHERE  THE  FIRST  SUBSCRIPT  REFERS  TO  THE 

C  HYPERCUBE  COORDINATES  IN  THE  SECOND  SUBSCRIPTS 

C  DIRECTION. 

C  VAL  -  THE  FUNCTIONAL  VALUES  AT  THE  CORNERS  OF  THE  HYPERCUBE 

C  SURROUNDING  XI.  THIS  VECTOR  MUST  BE  FILLED  EQUIVALENT 

C  TO  AN  NDIM  ARRAY  WITH  EACH  DIMENSION  AS  2.  THE  SIZE 

C  OF  VAL  SHOULD  BE  ATLEAST  2  * *N D I M . 

C  VM IN  -  A  MINIMUM  VALUE  OF  VAL  FOR  WHICH  THE  INTERPOLATION 
C  WILL  USE  A  CORNER  VALUE. 

C  WORK  -  A  WORK  VECTOR  OF  ATLEAST  NDIM*2.  USE  TO  STORE  CO0R- 
C  DINATE  WEIGHTS. 

C 

c 

C  OUTPUT 

C 

C  RETURNS  INTERPOLATED  VALUE  OF  VAL  AT  XI 

C 

C  CALLED  BY  MOMENT 

C 

C  *************************************************************** 

D I  MENS  ION  XI  (1  ) , XVAL ( 1  ) , VAL ( 1  ) , WORK ( 1  ) 

C 

C  SET  UP  THE  COORDINATE  WEIGHTS 
C 

NDI-I ABS (NDIM) 

IF (NDIM  . LT.  0)  GO  TO  1 
DO  100  I-l.NDI 
I  2  "I  *2 
I  1 -I  2-1 

WORK  (1  2  )  -  (XI  ( I  )  -  XVAL  ( I  1  )  )  /  (X  VAL  ( I  2  ) -XVAL  (  I  I  )  ) 

WOR  K ( I  1  )-l.  -  W  OR  K ( I  2) 

100  CONTINUE 

INTERPOLATE  -  USE  BINARY  COUNTER  FOR  COORDINATE  LOCATION 
C 

1  DTERPI-0. 

SUM-0. 


'  n 


i 
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000.  oooooooooo 


ND-2**NDI 
DO  201  I-l.ND 

IF ( VAL ( I )  .LT.  VM1N)  GO  TO  201 
L  -1  -  1 
WEIGHT-  1. 

DO  200  J- 1 , NDI 
N -M OD (L , 2  )  +  J  *2  -  1 
WE1CHT-WEIGHT*W0RK(N ) 

L-L  /2 

200  CONTINUE 

SUM -SUM  +  WEIGHT 

DTERP  I-DTERPI  +  W EIGHT *VAL ( I  ) 

201  CONTINUE 

IF  (SUM  .  EQ.  0.  )  GO  TO  202 
DTERP I-DTERP  I  /  SUM 
RETURN 
C 

2  C  2  STOP 
END 


C 

C 

c 


*  *  1 


SUBFILE  2  8 


SUBROUTINE  GREEN (Z , Z 1 , T .ALPHA, TO, IER> 

IMPLICIT  INTEGER**  (I-NJ 

********************************************************** 


PUR  RISE 

TO  COMTUTE  THE  GENERALIZED  CREENS  FUNCTION 
USES  GREEN  1 

SEE  GREEN  1  FOR  ARGUMENT  LIST 


************************************************************ 

REAL  N,N 

COMMON  /WNDPRM/D XZ 0 , DY XO , DZO , U 0, M , N , Z I K V 

IF  (I.  .Eg.  1.  )  GO  TO  2 

X  ^  =  n  ' N 

AT-ALPHA*! 

TC-O. 

IF (AT  .CE.  Z1 )  RETURN 

CALL  CF.EEN1  (  (Z-hAT  )  **X  2,  Z  J  **X2,  X2*X2*T ,  (N- 1 .  )  /X  2,  71 ,  I ER  ) 

T 1-7  1 *X  2*Z 1 ** ( 1 . -N ) 

U«  1 . 

T2-0. 

IF (ABS (ALPHA )  .LT.  l.E-4)  GO  TO  1 
ZMZ-Z-Z 1 +AT 
X2-N  +  1. 

A!.1-ALP11A*X2 

ZMZN-Z  1**X2  -  (Z 1  -AT ) **X2 
ARG-  (-AN  1 *ZMZ*ZMZ )  /  (4. *ZMZN) 

IF  (ARC.  .LT.  -70.  )  GO  TO  3 
T 2 -SORT (AN  1/ (4.  *3.  1415926*ZMZN) )*EXP(ARC) 


nnnnnnnnnnnnnnnnnnnnnnnonnnnnnnnnnnnnn  nnn  non 


3  IF  (T  I.  LT.  1.  E-30  .AND.  T2.  LT.  1.  E-30)  RETURN 

CALCULATION  OF  MIXING  RATIO,  U,  BY  N-l  ANALOCY 

CALL  CREEN1 (Z4AT ,Z1,T,0. ,T1U,IER) 

X2=2. 

AN  1  «=ALPHA*X  2 

ZMZN=Z1**X2  -  (Z1-AT)**X2 
T2H-0. 

ARC  « ( -AN  1 *ZMZ*  ZMZ ) / ( 4 . *  ZMZN ) 

IF (ARG  .LT.  -70. )  GO  TO  4 

T2U*=SQRT (AN  1/(4.  *3.  1415926*ZKZN) )*EXP(ARC) 

4  IF (T 1U. LT. 1.  E-30  .AND.  T2U.  LT.  1 .  E-30)  GO  TO  1 
CALL  GREEN 1 (Z,Z1,T,ALFHA,G, 1ER) 

U  *=  (G  -T  2U  )  /  (T  1U  -T  2U ) 

1  I F (U  .LT.  0.  )  U>=0. 

IF (U  .CT.  1.  )  U«l. 

COMBINE  LIMITING  SOLUTIONS  WITH  DETERMINED  MIXING  RATIO 

T0HJ*T  1  +  ( ! .  -U ) *T 2 
RETURN 

2  CALL  CREEN1  (Z,Z1,T, ALPHA, TO,  IER) 

RETURN 

END 

*********************  SUBFILE  29  ************************** 


SUBROUTINE  GREEN  1 
PURPOSE 

COMPUTE  THE  I  BESSEL  FUNCTION  FOR  A  CIVEN  ARGUMENT  AND  ORDER 
AND  MULTIPLY  BY  AN  APPROPRIATE  POWER  OF  THE  ARCUMENT 
AND  AN  EXPONENTIAL  IN  ORDER  TO  CALCUUTE  THE  CREENS 
FUNCTION  FOR  THE  WIND  DIFFUSION  EQUATION 


USAGE 

CALL  GREEN1 (Z,Z1,T,NU,EI, IER) 

DESCRIPTION  OF  PARAMETERS 

Z , Z 1 , T  -THE  ARGUMENTS  OF  THE  FUNCTION  DESIRED 
HU  -THE  ORDER  OF  THE  I  BESSEL  FUNCTION 
BI  -THE  RESULTANT  BESSEL  FUNCTION 
IER  -RESULTANT  ERROR  CODE  WHERE 

IER— 1  EXPONENTIAL  UNDERFLOW  (NON-FATAL),  BI  SET  TO  0.  C 
IEK-0  NO  ERROR 

IER-1  NU  NEAP  NEGATIVE  INTEGER 
IER*=2  OVERFLOW'  IN  GA MMA 

1ER-3  UNDERFLOW,  BI  .LT.  l.E-32,  BI  SET  TO  0.0 
I ER“4  OVERFLOW,  X  .CT.  90  WHERE  X  .GT.  N 
IER-5  X  IS  NEGATIVE 

REMARKS 

NU  ISA  REAL  NUMBER 
N  AND  X  MUST  BE  .GE.  ZERO 

THIS  SUBROUTINE  IS  A  MODIFICATION  OF  BESI  WHICH  COMPUTES  THE 
I  BESSEL  FUNCTION  FOR  INTEGER  ORDERS.  THE  CHANCE  REQUIRES  THE 
USE  OF  THE  GAMMA  FUNCTION  FOR  COMPUTING  THE  FIRST  TERM  OF  THE 
SERIES.  THE  SUCCESSIVE  TERMS  ARE  CALCULATED  WITH  THE  SAME 
RECURSION  FORMULA  AND  THE  ASYMPTOTIC  APPROXIMATION  IS  ALSO 


i  < 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


UNCHANGED.  BESI  IS  IN  THE  IBM  SYSTEM/360  SCIENTIFIC  SUBROUTINE 
PACKAGE.  MODIFICATIONS  MADE  BY  D.  DWRE,  AERODYNE  RESEARCH, 
INC.  JANUARY  15,  1979. 

SUBROUTINES  AND  FUNCTIONS  REQUIRED 

GAMMA  WHICH  COMPUTES  THE  GAMMA  FUNCTION 

METHOD 

COMPUTES  I  BESSEL  FUNCTION  USING  SERIES  OR  ASYMPTOTIC 
APPROXIMATION  DEPENDING  ON  THE  RANGE  OF  THE  ARGUMENT. 

CALLED  BY  MOMENT 


SUBROUTINE  GREEN  1 (Z , 2 1 , T, NU, BI , I ER ) 

IMPLICIT  INTEGERS  (I-N) 

REAL  NU 

X-2.  *S QRT  (Z*Z  1 )  /T 

CHECK  FOR  ERRORS  IN  NU  AND  X  AND  EXIT  IF  ANY  ARE  PRESEN! 

IER-0 
B I*  1 . 0 

IF (NU ) 1 0, 15,10 
10  IF(X ) 1 60, 20, 20 
15  I F (X ) 1 60, 17,20 
17  ARG-- (Z+2  1 )  /T 

IF (ARC  .LT.  -80. )  GO  TO  170 

B  I-EXP  (ARC  )  /T 

RETURN 

DEFINE  TOLERANCE 
20  TOL-l.E-3 


IF  ARGUMENT  GT  12  AND  CT  NU,  USE  ASYMPTOTIC  FORM 

IF  (X -1 2.  ) 4 0 ,  40.  ?r 
30  1  F(X  -ABS  (NU  )  )4C,  40,  1  1 C 

COMPUTE  FIRST  TERM  OF  SERIES  AND  SET  INITIAL  VALUE  OF  THE  SUM 


40  XX-X/2. 

N-IN7 (NU  ) 

FN-N 

R-NU-FN 

CALL  GAMMA(1.4NU,CR,IER) 
IF( IER  .EQ.  0)  GO  TO  60 
50  BI-0.0 
RETURN 

60  TERM-l./GR 
70  B I -TERM 
XX-XX*XX 


COMPUTE  TERMS,  STOPPING  WHEN  ABS (TERM)  LE  AES (SUM  OF  TERMS ) *TOLERANC E 
DO  90  K-l, 1000 

IF (ABS (TERM) -ABS(B1*T0L))95, 95, 80 
60  FK-K 

FK-FK*  (NU+FK ) 

TERM-TERM* (XX/FK) 

90  BI-BI+TERM 
95  ARC—  (Z+Z1)/T 
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c 

c 

c 


IF (ARC  .LT.  -80.)  CO  TO  170 
BI«BI*(Z1/T) **NU*EXP(ARC)  /I 

RETURN  B I  AS  ANSWER 


100  RETURN 
C 

C  X  CT  12  AND  X  CT  NU,  SO  USE  ASYMPTOTIC  APPROXIMATION 
C 

110  FN-4.*NU*NU 
115  XX-1./(8.*X) 

TERM*  1. 

BI*1. 

DO  130  K=l, 30 

IF(ABS  (TERM)  -ABS  (BI*TOL)  )  UO,  140,  120 
120  FK=  (2*K-1 )  **2 

T  ERM-TERM*XX*(FK-FN)  /FLOAT  (K  ) 

130  B I =B 1 +T  ERM 
C 

C  SIGNIFICANCE  LOST  AFTER  30  TERMS,  TRY  SERIES 
C 

CO  TO  40 

140  Pl-3. 1 4 1 5S2653 
ARC=X- (Z+Z 1 ) /T 
IF (ARC  .LT.  -SO. )  GO  TO  170 

B 1*B I* (Z 1 /Z) ** (NU/2.  )*EXP(AKC)/SQRT(2.*PI*X) /T 
GO  TO  100 
160  IER  =  5 

GO  TO  100 
170  B  1*0. C 
GO  TO  50 
END 
C 

C  ******************  SUBFILE  30  **************************** 

C 

C 

C  . 

c 

C  SUBROUTINE  GAMMA 

C 

C  PURPOSE 

C  COMPUTES  THE  GAMMA  FUNCTION  FOR  A  CIVEN  ARGUMENT 

C 

C  USAGE 

C  CALL  GAMMA.  (XX,  CX,  I ER  ) 

C 

C  DESCRIPTION  OF  PARAMETERS 

C  XX  -THE  ARGUMENT  FOR  THE  GAMMA  FUNCTION 

C  GX  -THE  RESULTANT  GAMMA  FUNCTION 

C  I ER  -THE  RESULTANT  ERROR  CODE  WHERE 

C  I  ER=0  NO  ERROR 

C  I EP.  =  1  XX  IS  WITHIN  .000001  OF  BEING  A  NEGATIVE  INTEGER 

C  I  ER  =  2  XX  CT  57,  OVERFLOW,  CX  SET  TO  1.E32 

C 

C  COMMENTS 

C  NONE 

C 

C  SUBROUTINES  AND  FUNCTIONS 

C  NONE 

C 

C  METHOD 

C  THE  RECURSION  RLATION  AND  POLYNOMIAL  APPROXIMATION 

C  BY  C.  HASTINGS, JR. .  'APPROXIMATIONS  FOR  DIGITAL  COMPUTERS', 

C  PRINCETON  UNIVERSITY  PRESS,  1955 


<  • 


t 


<-!  on  o  o  nn  on  n  n  r:  o  no  C!<~!  nr-  r~  ^  o  r-  n  nno  non 


SUBROUTINE  CAMMA(XX,  GX,  IER  ) 

IMPLICIT  INTEGER*^  (I-N) 

IF (XX-5  7.  )  6,6,4 
4  I  ER  =  2 
GX-1.E32 
RETURN 
6  X*XX 

ERR*1. OE-6 
I  EK  =  0 
CX=1.  C 

IF(X-2. 0)  50,  50,  1  5 
10  IF (X-2. 0)  1  10, 110,15 
15  X-X-1.0 
C  X=G  X  *X 
GO  TO  10 

50  IF(X-l.C)  60, 120, 110 

SEE  IF  X  IS  NEAR  NEGATIVE  INTEGER  OR  ZERO 

60  IF(X-ERR)  62,62,80 
62  Y  =F LOAT (I  NT (X  ) )  -X 

I F  (ABS  (V  )  -ERR. )  1  30,  130,  70 

X  NOT  NEAR  A  NEGATIVE  INTECER  OR  ZERO 

7C  I  F ( X - 1 . 0)80,  8.0,  110 
K  CXCX/X 
X-X  +  l.  0 
GO  TO  7C 
11C  Y=X-1. 0 

GY«1.0+Y*(-0.  5771C1  7+Y*(0.  9658540+Y*  (-0.  8  764  2 1  8+Y  *  (0.  8  3282 1  2+ 
1 Y* (-0.  5684729+Y*(0.  254  8205+Y* (-0.  05 1 49930) )))))) 

G  X  *G  X  *G  5 
1  2(  RETURN 
1 2u  IER*  1 
RETURN 
END 

S  UB ROUT  I N  E  RKM  (N  , XL , XU  ,  Y,  HM  IN ,  DEL, ACCUR C , WK, ND ) 

IMPLICIT  INTECER*4  (I-N) 


NUMERICAL  INTEGRATION  ROUTINE  FOR  SYSTEMS  OF  ODE'S 
USING  THE  RUNCE-KUTTA-MERSON  TECHNIQUE 


INPUT  PARAMETERS 


N 

XL 

XU 

Y 


HMIN 

DEL 


C. 


NUMBER  OF  FIRST  ORDER  DIFFERENTIAL  EQUATIONS 
INITIAL  ABCISSA  OF  TLI  INTERVAL 
THE  FINAL  ABCISSA  OF  THE  INTEGRATION  INTERVAL 
A  SINGLY  DIMENSIONED  ARRAY  OF  LENGTH  N.  WHEN 
RKM  IS  CALLED  IT  MUST  CONTAIN  THE  VALUES  OF 
THE  DEPENDENT  VARIABLES  AT  XL.  UPON  RETURN 
TO  THE  CALLING  PROGRAM  Y  CONTAINS  THF  VALUES 
OF  THT  DEPENDENT  VARIABLES  AT  XT. 

THE  MINIMUM  STEP  SIZE  THAT  WILL  BE  USED  FOR  THE 
INTEGRATION 

THE  INITIAL  ESTIMATE  OF  THE  STEP  SIZE  AND  UPON 
RETURN  TO  THE  CALLING  PROGRAM  DEL  CONTAINS  THE 
FINAL  STEP  SLZE  USED.  THIS  VALUE  SHOULD  BE  USED 
IN  THE  NEXT  CALL  TO  PRODUCE  AN  EFFICIENT  INTEGRATION. 
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C  DEL  IS  RETURNED  WITH  THE  VALUE  ZERO  IF  IT  HAS 

C  BEEN  HALVED  BELOW  HMIN. 

C  ACCURC  -  PREASSIGNED  ACCURACY  WHICH  IS  ALSO  USED  IN  ADJUSTING 
C  THE  STEP  SIZE. 

C  WK  -  AT  LEAST  A  BLOCK  OF  N  BY  6  FLOATING  POINT  LOCATIONS 
C  USED  FOR  A  WORK  ARRAY. 

C  ND  -  THE  DIMENSION  OF  ARRAYS  Y  AND  WK. 

C 

C 

C  IT  IS  REQUIRED  THAT  THE  USER  OF  RKM  WRITE  A  SUBROUTINE 
C  DEFINING  THE  DIFFERENTIAL  EQUATIONS.  THE  SUBROUTINE 
C  STATEMENT  SHOULD  LOOK  LIKE  -  SUBROUTINE  DIFEQ (N ,X, Y, YP)  . 

C 

C  WHERE 

C  N  -  THE  NUMBER  OF  EQUATIONS 

C  X  -  THE  INDEPENDENT  VARIABLE 

C  Y  -  SINGLY  DIMENSIONED  ARRAY  OF  DEPENDENT  VARIABLES 

C  YP  -  SINGLY  DIMENSIONED  ARRAY  OF  THE  RATES  CF  Y  AT  X 

C  YP(I  )  -  D  Y (I  )/DX 

C 
C 

DIMENSION  Y  (ND  )  ,WK(t.'D  ,  6  ) 

LOGICAL  FIRST, QUIT 
C 

C  SET  UP  NEEDED  VARIABLES  UPON  ENTRY 
C 

XN-XL 

H=DEL 

FIRST-. TRUE. 

QUIT-. FALSE. 

C 

C  CHECK  IF  XN  IS  CLOSE  TO  XU 
C 

20  I F  (XN-H!  .LT.  XU)  CO  TO  30 
DEL-1! 

H-XU-XN 
QUIT-. TRUE. 

IF ( FIRST  )  DEL-ii 
C 

C  MAKE  FIRST  CALL  TO  DIFEQ  AT  THE  BEGINNING  OF  INTERVAL 
C 

30  CALL  DIF EQ (N ,XN , Y, VK ( 1 , 1 ) ) 

C 

C  PERFORM  THE  RUNGE-KUTTA-MERSON  ALGORITHM 
C 

40  H3^!/3. 

DO  5C  I  -  1 , N 

WK (I , 3 )=H  3*WK (1,1) 

50  WK ( I  ,  6  )-Y  (I  )  +W K  ( I  ,  3 ) 

CALL  DlF EQ  (N  ,  XN-Hi  j,  WK  ( 1 , 6 ) ,  WK  ( 1 , 2 )  ) 

DO  60  1-1, N 

60  WK ( 1 , 6 )=Y (I )+  (WK (I , 3 )+H 3*WK(1 ,2) )/2. 

CALL  D1 F  EQ  (N  ,  XN-HI  3,  WK(  1 , 6  )  ,WK  ( 1 , 2  )  ) 

DO  70  1-1/ N 

WK ( I ,  4 )  -H 3*WK ( I ,2) 

70  WK (1 , 6 )«Y ( I )  +  (3. *WK ( 1 , 3)+9. *WK(I ,4 ) ) /8. 

CALL  DIFEQIN  .XN-H1/2.  ,  WK  ( 1 , 6  ) ,  WK  ( 1  ,  2  )  ) 

DO  EO  I = 1 , N 

WK ( I , 5 )«H 3*WK(I ,2) 

80  WK(I  ,  6 )-Y  (1  )+  (3.  *WK (1,3 ) -9.  *W'K (I  ,  4  )  + 1  2.  *WK (I  ,  5  )  )  /2 . 

CALL  DI FEQ (N ,XN+H , WK ( 1 , 6 ) , WK ( 1 , 2 ) ) 

C 

C  FIND  THE  LARGEST  RELATIVE  ERROR 
C 
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TLST-O. 

DO  9 j  I-l.N 


VX=Y  (1  ) 

IF  ( Y X  .  EQ. 
E= ( (Wl 
90  TEST 


n  ^  y  vaArrimr 
K(K3)-9.*WK(I  ,«)/2.+4.*WK(I,5)-H3*WK(I  ,2)/2.  )  fb. 
AMAa 1 (TEST , ABS (E ) ) 


)  /YX 


FIRST". FALSE. 

IF (TEST  .LT.  AC CUR C)  GO  TO 


100 


IF  THE  LARGEST  ERROR  IS  GREATER  THAN  ACCL'RC  HALF  THE  STEP 
SIZE  AND  TRY  AGAIN. 


H  =H  /  2 . 

I F (H  .LT.  HMIN)  GO  TO  10 
QUIT". FALSE. 

CO  TO  A l 

TRUNCATION  ERROR  LESS  THAN  ACCL’RC,  RES  FT  THF  Y  ARRAY  TO 
C  SET  UP  FOR  THE  NEXT  INTERVAL 

100  XN=XN-K 

DO  110  I  =  1 ,  N 

lit  Y  (1  )=Y  (1  )  +  (WK  (I  ,  3 )  +A .  *WK  ( 1 , 5  )+H  3*VK  (I,2))/2. 

C 

C  CHECK  FOR  STEP  SIZE  DOUBLING.  DOUBLE  IF  LARGEST  REL/  *TVE 
C  ERROR  IS  32  TIMES  LESS  THAN  ACCL'RC. 

IF <  .NOT. (TEST  .CE.  ACCLRC/32.  .OR.  QUIT))  H"H+H 
IF  (.NOT.  QUIT)  Gt.  TO  2( 

RET  UK N 

c 

C  THE  VALUE  OF  H  (DEL)  IS  LESS  THAN  THE  SPECIFIED  MINIMUM. 

C  RF.KF.I  THIS  AND  ERROR  OIT. 

C 

10  CONTINUE 
Dl,L=t  . 

RETURN 

LN. 

ROT  Tor 
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